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Abstract. Supernova models with a full spectral treatment of the neutrino transport are presented, employing the 
Prometheus/Vertex neutrino-hydrodynamics code with a variable Eddington factor closure of the 0{v/c) moments equations 
of neutrino number, energy, and momentum. Our "ray-by-ray plus" approximation developed for two- (or three-) dimensional 
problems assumes that the local neutrino distribution function is azimuthally symmetric around the radial direction, which 
implies that the nonradial flux components disappear. Other terms containing the angular velocity components are retained 
in the moments equations and establish a coupling of the transport at different latitudes by lateral derivatives. Also lateral 
components of the neutrino pressure gradients are included in the hydrodynamics equations. This approximative approach for 
neutrino transport in multi-dimensional environments is motivated and critically assessed with respect to its capabilities, limi- 
tations, and inaccuracies in the context of supernova simulations. In this first paper of a series, one- (ID) and two-dimensional 
(2D) core-collapse calculations for a (nonrotating) 15 Mq star are discussed, uncertainties in the treatment of the equation 
of state - numerical and physical - are tested, Newtonian results are compared with simulations using a general relativistic 
potential, bremsstrahlung and interactions of neutrinos of different flavors are investigated, and the standard approximation in 
neutrino-nucleon interactions with zero energy transfer is replaced by rates that include corrections due to nucleon recoil, ther- 
mal motions, weak magnetism, and nucleon correlations. Models with the full implementation of the "ray-by-ray plus" spectral 
transport were found not to explode, neither in spherical symmetry nor in 2D when the computational grid is constrained to a lat- 
eral wedge (< +45°) around the equator. The success of previous two-dimensional simulations with grey, flux-limited neutrino 
diffusion can therefore not be confirmed. An explosion is obtained in 2D for the considered 15 Mq progenitor, when the radial 
velocity terms in the neutrino momentum equation are omitted. This manipulation increases the neutrino energy density in the 
convective gain layer by about 20-30% and thus the integral neutrino energy deposition in this region by about a factor of two 
compared to the non-exploding 2D model with the full transport. The spectral treatment of the transport and detailed description 
of charged-current processes leads to proton-rich neutrino-heated ejecta, removing the problem that previous explosion models 
with approximate neutrino treatment overproduced N = 50 closed neutron shell nuclei by large factors. 
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1. Introduction 

Convective processes and hydrodynamic instabilities play an 
important role in different regions of collapsing and explod- 
ing stars. (1) Inside the nascent neutron star, i.e. below and 
around the neutrinosphere, convection enhances the transport 
and release of neutrinos and thus can raise the neutrino lu- 
minosities with potentially helpful consequences for the de- 
layed shock revival and the neutrino-driven explosion mecha- 
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nism (e.g.. Burrows 1987; Burrows & Lattimer 1988; Miralles 
et al. 2000; Pons et al. 1999; Wilson & Mayle 1988, 1993; 
Keil et al. 1996; Keil 1997). (2) In the neutrino-heating region 
between gain radius and stalled supernova shock, convective 
overturn was found to enhance the energy deposition by neu- 
trinos and the energy transport to the shock. Therefore it al- 
lows for neutrino-driven explosions even if models in spherical 
symmetry fail (Herant et al. 1994; Burrows et al. 1995; Janka 
& Miiller 1996). (3) Multi-dimensional processes in the outer 
layers of the exploding star destroy the onion-shell structure 
of the progenitor and are responsible for the anisotropics and 
macroscopic mixing of chemical elements which was observed 
in Supernova 1987 A and other supernovae, and which possi- 
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bly is the origin of many morphological features of supernova 
remnants (e.g., Kifonidis et al. 2003 and references therein). 

The existence of convectively unstable stratifications in the 
supernova core had long been suspected on grounds of spheri- 
cally symmetric models (e.g., Epstein 1979; Bethe et al. 1987; 
Burrows 1987; Bethe 1990) before confirmation became possi- 
ble by multi-dimensional modeUng with sufficient resolution. 

The first such multi-dimensional simulations (Burrows & 
Fryxell 1992; Janka & Mtiller 1993; Miiller 1993) still ignored 
neutrino efl'ects, were able to follow the evolution for a very 
limited period of time only, and started from post-bounce mod- 
els that were constructed artificially or were adopted from sim- 
ulations based on different assumptions about the input physics. 
Therefore they were only suitable to demonstrate the possi- 
bility in principle, but otherwise showed the transient decay 
of an initially unstable postshock stratification and thus could 
not be considered as meaningful models of supernovae (cf. 
Bruenn & Mezzacappa 1994). However, simulations including 
the effects of neutrino heating and cooling (Herant et al. 1992) 
and true core collapse environments (Janka & Miiller 1994, 
1995, 1996), neutrino transport by grey, flux-limited diffusion 
(Herant et al. 1994; Burrows et al. 1995), multi-group ID trans- 
port (Mezzacappa et al. 1998b), rotation (Fryer & Heger 2000), 
and three-dimensional hydrodynamics (Fryer & Warren 2002, 
2004) confirmed robustly the development of convective over- 
turn in the neutrino-heated layer just behind the stalled shock, 
although it is still a matter of debate whether the associated 
convective effects are sufficient for driving delayed explosions 
(Mezzacappa et al. 1998b). 

The situation is more controversial for convective activity 
below the neutrinosphere. A detailed investigation of neutrino 
transport properties at the conditions inside nascent neutron 
stars (Bruenn et al. 2004) revealed no evidence for neutron 
finger instabiUty as assumed by the Livermore group (Wilson 
& Mayle 1988, 1993) to trigger neutrino-driven explosions by 
convectively enhanced neutrino fluxes. In the simulations of 
Janka & Miiller (1996) the use of an inner boundary condi- 
tion with imposed neutrino luminosities certainly played a role, 
causing neutrino absorption and heating and thus convective 
energy transport in the neutron star layers around the inner 
boundary (neutrino diffusion was ignored!). This convection 
below the neutrinosphere extended outward to a region where 
cooling by neutrino losses was fast enough to take up the en- 
ergy delivered by the convective transport. However, full-scale 
hydrodynamic models of nascent neutron stars in two dimen- 
sions including a simple description of grey, radial neutrino dif- 
fusion and its effects on lepton number and energy transport 
showed the rapid development of convection inside the proto- 
neutron star (Burrows & Fryxell 1993) and its persistence in 
a growing volume until at least one second in nonrotating as 
well as rotating cases (Keil et al. 1996; Keil 1997; Janka & 
Keil 1998; Janka 2004, whose simulations covered the evolu- 
tion until at most 1.2 seconds after bounce). A careful analy- 
sis (see Keil 1997) revealed that the exchange of energy and in 
particular of lepton number between moving fluid elements and 
their surroundings was important so that the convective mode 
could not be described by ideal Ledoux convection. Instead, it 
had the character of a doubly diflxisive instabiUty, which can 



occur at Ledoux stable conditions as discussed by Bruerm & 
Dineva (1996); Bruenn et al. (1995) and more recently in much 
detail by Bruenn et al. (2004). 

This importance of entropy and lepton exchange of buoyant 
fluid elements with their surroundings via neutrino transport 
was also pointed out by Mezzacappa et al. (1998a). The latter 
authors, moreover, argued on the basis of their numerical and 
analytic models and timescale analyses that fast equilibration 
by neutrino transport tends to damp or even suppress entropy- 
and lepton-driven proto-neutron star convection around and 
somewhat below the neutrinosphere, making it inefficient in the 
transport of lepton number and energy. The authors, however, 
admit that their investigations lack multi-dimensional "realis- 
tic" multi-group neutrino transport, coupled self-consistently 
to multi-dimensional hydrodynamics, which have to be applied 
in time-dependent evolution simulations of supernova cores be- 
fore final conclusions can be drawn about the development of 
convection in optically thick regions. Bruenn et al. (2004) in 
particular stress the importance of a careful inclusion of neu- 
trino transport effects for a reliable description of doubly diffu- 
sive instabihties. 

While the basic results of the models of Keil et al. (1996) 
might therefore still withstand critical assessment by future, 
better simulations, quantitative aspects (e.g., the speed at which 
the convective layer grows inward towards the stellar center and 
enhances the deleptonization and energy loss) were undoubt- 
edly affected by the simplicity of the treatment of the neutrino 
physics. Keil et al. (1996) used the assumption that neutrinos 
are in equilibrium with the stellar fluid (neutrino energy and 
pressure were therefore added to the corresponding gas quanti- 
ties), that neutrino transport can be approximated by first-order 
grey, (flux-limited) radial equilibrium diffusion (cf. Mihalas & 
Mihalas 1984, §97), and that the diffusive ffuxes of lepton num- 
ber and energy can be combined for electron neutrinos and an- 
tineutrinos into a gradient form with a single, average diffusion 
coefficient (cf. Appendix A of Keil & Janka 1995). The need 
for removing these deficiencies has motivated new efforts for 
improved modeling, which were recently reported by Swesty 
& Myra (2005). 

In addition to Ledoux convection and doubly diffusive in- 
stabilities another form of multi-dimensional hydrodynamic in- 
stability seems to play a role in the supernova core: While 
one-dimensional, time-dependent hydrodynamic simulations 
(without neutrino transport) found that the standing accretion 
shock is stable against radial perturbations, two- and three- 
dimensional (adiabatic) simulations demonstrated that it is un- 
stable to non-radial modes, which lead to rapid growth of tur- 
bulence behind the shock and thus shock expansion and de- 
formation (Blondin et al. 2003; Mezzacappa & Blondin 2003). 
Such a standing accretion shock instability, SASI, was actu- 
ally predicted for accreting black holes by analytic considera- 
tions by Foglizzo (2001, 2002), who described its origin from 
an amplifying feedback cycle of entropy and vorticity pertur- 
bations, which are advected inward and trigger acoustic waves 
that propagate outward and distort the shock, thus closing the 
cycle. The existence of this phenomenon at non-adiabatic con- 
ditions in the postshock flow was confirmed by supernova sim- 
ulations including a simplified treatment of neutrino transport 



R. Buras et al.: Two-dimensional hydrodynamic core-collapse supernova simulations witii specb-al neutiino ti'ansport 



3 



(Scheck et al. 2004; Scheck 2005). Its action can clearly be 
seen in particular if neutrino heating is too weak to drive strong 
postshock convection. In agreement with the adiabatic mod- 
els of Foglizzo (2001, 2002) and Blondin et al. (2003), Scheck 
et al. (2004) find highest growth rates for low (I - 1,2) modes, 
which produce large, global asymmetries of the supernova ex- 
plosion and might explain the observed pulsar kick velocities 
and grossly asymmetric distribution of elements seen in super- 
nova remnants Uke Cassiopeia A (see Janka et al. 2005c). 

Advancing core collapse supernova modeling to the next 
level of complexity and reliability beyond the stage of merely 
demonstrating possibiUties in principle, requires improve- 
ments in the handling of neutrino transport coupled to multi- 
dimensional hydrodynamic simulations. Such improvements 
are necessary to clarify which role hydrodynamic instabilities 
play in the supernova core and how these are linked to the 
mechanism of the explosion. Crucial questions in this context 
are: Does convection take place below the neutrinosphere and 
if so, which kind of instability is it? Does proto-neutron star 
convection raise the neutrinospheric luminosities and thus the 
neutrino heating behind the shock on timescales relevant for 
the explosion? Does neutrino-driven convection in the heat- 
ing layer become sufficiently strong to trigger the explosion? 
Can the positive results of previous simulations with simpUfied 
neutrino transport by grey, flux-limited difliision (Herant et al. 
1994; Fryer 1999; Fryer & Heger 2000; Fryer & Warren 2002, 
2004) be confirmed? Is the accretion shock instability found by 
Blondin et al. (2003) important and how does it interact with 
convective instabilities and rotation? 

With the aim to address these questions in a next genera- 
tion of core collapse studies, we have developed a new method 
for treating time- and frequency-dependent neutrino transport 
in multi-dimensional (in this work two-dimensional) environ- 
ments. Our transport code is a generalization of the ID ver- 
sion Vertex of Rampp & Janka (2002), adopting the approx- 
imations to general relativity described in the latter paper and 
tested for spherically symmetric problems against fully rela- 
tivistic calculations in Liebendorfer et al. (2005) and Marek 
et al. (2005). The code employs a detailed spectral description 
of neutrino-matter interactions with fully implicit energy-bin 
coupling and with a set of processes which was updated from 
Rampp & Janka (2002) by including interactions between neu- 
trinos of different flavors and improving neutrino-nucleon in- 
teractions such that the effects of nucleon thermal motions, re- 
coil, weak magnetism, and correlations in charged and neutral 
current processes are taken into account. The high dimension- 
ality of the spatial and momentum dependence of the phase 
space distribution function is reduced by assuming the neu- 
trino intensity to be axially symmetric around the radial direc- 
tion. Consistent with this, we ignore lateral flux components 
and effects of neutrino viscosity and come up with a set of mo- 
ments equations for neutrino number density, energy density 
and radial flux, which include all velocity-dependent terms of 
0{vjc) for radial and lateral velocity components. These equa- 
tions are solved in each angular bin of the spatial grid, tak- 
ing into account the lateral (and in 3D also azimuthal) cou- 
pling by corresponding derivative terms in an operator splitting 
step ("ray-by-ray plus" approach). Closure of the set of mo- 



ments equations is achieved by using variable Eddington fac- 
tors, which are computed from a model Boltzmann equation 
as described in Rampp & Janka (2002). The efliciency of the 
code performance is increased by determining these normal- 
ized moments of the neutrino intensity, which usually exhibit 
relatively weak variability with lateral direction, only once for 
all latitudes on grounds of an angularly averaged stellar back- 
ground. Because of the fact that we refer to the solution of 
the model Boltzmann equation for closing the set of moments 
equations, we named our code MuDBaTH (Multi-Dimensional 
Boltzmann Transport and Hydrodynamics). 

Our approach to deal with neutrino transport in two- (and 
in an analogous way in three-) dimensional environments is in- 
tended to be a valid approximation for situations where only 
local macroscopic inhomogeneities are present, e.g. due to con- 
vection. It is, however, likely to become poor when a large 
global deformation (e.g. by rapid rotation) develops and the 
neutrino flow is therefore not well described by an at least 
on average dominant radial component. Our approach is com- 
plementary concerning the employed approximations to that 
of Livne et al. (2004), who introduced a scheme for solving 
the frequency-dependent Boltzmann transport equation in ax- 
ially symmetric problems, but sacrificed energy-bin coupling 
and the majority of the terms with dependence on the motion 
of the stellar plasma. It also differs from the two-dimensional 
multi-group flux-limited diffusion approximation employed by 
Walder et al. (2005), who abandon energy-bin coupling and 
velocity-dependent terms, too. 

In this Paper I and a subsequent one (Buras et al. 2005; 
Paper II) we present the first two-dimensional hydrodynamic 
simulations with a full spectral implementation of neutrino 
transport, applied to the core collapse and post-bounce evolu- 
tion of progenitor stars with different masses (in the range from 
11.2 Mq to 25 Mq), compared to corresponding simulations in 
spherical symmetry. In Paper I we shall describe in detail our 
method for neutrino transport in two spatial dimensions, dis- 
cuss numerical aspects and tests of the code and input physics, 
and try to assess the limitations of our method. We will con- 
centrate on the presentation of numerical results for the evo- 
lution of a 15 Mq progenitor star (Model sl5s7b2 of Woosley 
& Weaver 1995) with varied input physics. For the 2D (axi- 
ally symmetric) simulations we confine ourselves to using a 
roughly 90 degree angular wedge around the equatorial plane 
and will focus our discussion on neutrino-driven convection in 
the hot bubble region and its effect on the supernova dynamics. 
In Paper II we shall present and compare ID and 2D simula- 
tions for different progenitors, using different angular wedges 
and resolution (including a model with a full 180° polar grid) 
and considering also core rotation. Paper II will contain a de- 
tailed analysis of the importance of convection below the neu- 
trinosphere relative to convection in the neutrino-heating layer, 
and of the role of low-mode instabilities, pre-coUapse pertur- 
bations in the stellar core, and rotation in the investigated mod- 
els. Since our spectral treatment of neutrino transport in the 
2D simulations, which cover evolution periods of up to 300 ms 
after core bounce, is very CPU-time consuming - one model 
requires (2 - 5) x 10'^ floating point operations, depending on 
the angular resolution and wedge size -, computer resources 
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available to us prevented us from performing simulations with 
a systematic variation of important but not finally determined 
degrees of freedom, e.g., of the magnitude and distribution of 
angular momentum in the progenitor core. We therefore are 
forced to constrain our discussion on some selective cases. 

The present paper is organized as follows. Section 2 (sup- 
plemented by Appendices A-C) contains a description of the 
employed equations, numerical algorithm, and tests of the nu- 
merical scheme and input physics, in particular of the equa- 
tion of state. Section 3 presents results for ID (Sect. 3.1) 
and 2D (Sect. 3.2) core collapse and post-bounce simula- 
tions of a 15 Mq progenitor with Newtonian gravity and with 
our approximative treatment of the effects of relativistic grav- 
ity, varied neutrino-matter interactions, and neutrino transport 
with and without taking into account the velocity-dependent 
terms in the neutrino momentum equation. Thus retreating 
from our most complete implementation of neutrino trans- 
port by neglecting these velocity-dependent terms, we get 
cases with interesting behaviour, which once more demonstrate 
the sensitivity of the supernova evolution to supposedly sec- 
ondary changes of neutrino transport: In the ID simulation 
(Sect. 3.1.4) we observe a large-amphtude oscillation of the 
neutrino-decoupling layer between neutrinosphere and accre- 
tion shock, which grows in amplitude until it leads to an explo- 
sion; this reminds us of the k mechanism of vibrational insta- 
bility of stars. In the corresponding 2D simulation (Sect. 3.2.2) 
strong post-shock convection develops and triggers a neutrino- 
powered supernova explosion, whose properties will be com- 
pared to observational constraints. In contrast, both the ID and 
2D models with the full transport implementation (Sects. 3.1.1 
and 3.2.1, respectively) do not develop explosions. In Sect. 4 
we give a summary of our main results and draw conclusions. 



2. Governing equations and numerical 
implementation 

In order to be able to exploit symmetries of the problem we 
work in a system of spherical coordinates with radius r, lati- 
tudinal angle j?^, and longitudinal angle ip. For all simulations 
performed so far we have assumed that azimuthal symmetry 
holds with respect to the polar axis, but a generalization of the 
presented method to three-dimensional situations is straightfor- 
ward. 

Like the one-dimensional version documented by Rampp 
& Janka (2002) the algorithm relies on an operator splitting ap- 
proach which means that the coupled system of evolution equa- 
tions is processed in two independent steps, a hydrodynamic 
step and a neutrino-transport/interaction step. In each timestep 
these two steps are solved subsequently. 



2.1. Hydrodynamics 

For an ideal (i.e. nonviscous) fluid characterized by the mass 
density p, the radial, lateral, and azimuthal components of the 
velocity vector (v^, v,y, v^), specific energy e - e + ^{v^ + v^ + 
v^) with e being the specific internal energy, and gas pressure 
p, the Eulerian, nonrelativistic equations of hydrodynamics in 



spherical coordinates and azimuthal symmetry read: 



dt dr 



r sin § d{} 
1 d 
sin ?? dd^ 



^l + vl dp 50) ^ 

- P + ^ = -P-^ + QMr , 

r or or 
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VrVft - coti? \dp _ p 



+ -^ = --^ + Qm&, (3) 
r r otr r otr 



+ P- 



VrV^ + VtfV^COtl^ _ ^ 



(4) 



1 d 



^P \'~dr^~^ 'm] ^'^^^ ^rQur + V,fQM,i , (5) 

where O denotes the gravitational potential of the fluid, and 
Qm = iQur' Qutf) and Qe are the neutrino source terms 
for momentum transfer and energy exchange, respectively. 
Eqs. (1-5) are closed by the equation of state (EoS) which 
yields the pressure p for given p, e, and composition. In the 
case of nuclear statistical equilibrium (NSE) a third indepen- 
dent variable, the electron fraction Y^, is sufficient to character- 
ize the composition. For this variable the evolution is computed 
according to a conservation equation 

|- (pi's) (r^pYc Vr) + ^ (sin &pY^ v&) = , 

(6) 

where the source term Qi^/(p/ni\,y) is the rate of change of 
the net electron fraction (i.e. the number fraction of electrons 
minus that of positrons) due to emission and absorption of 
electron-flavor neutrinos, and Wby is the baryon mass. In case 
the medium is not in NSE, an equation hke Eq. (6) has to be 
solved also for the abundance of each nucleus k, Yk = nk/riby, 
using 

^ (pYk) + -2^- [r^pYk Vr) + ^ (sin §pYk v^) = Rk , 

ot or^ ' r^mij air 

(7) 

where and «by are the number density of nucleus k and the 
baryon number density, respectively, and Rk = p6Yk/6t, where 
6Yk/5t is a source term that describes the rate of composition 
changes by nuclear reactions for species k. 

For the numerical integration of Eqs. (1-7) we employ the 
Newtonian finite-volume code PROMETHEUS (Fryxell et al. 
1989, 2000), which was supplemented by additional problem 
specific features (Keil 1997) and the improvements described 
in Kifonidis et al. (2003). PROMETHEUS is a direct Eulerian 
implementation of the Piecewise Parabolic Method (PPM) of 
Colella & Woodward (1984). As a time-explicit, third-order in 
space, second order in time Godunov scheme with a Riemann 



R. Buras et al.: Two-dimensional hydrodynamic core-collapse supernova simulations witii specb-al neutiino ti'ansport 



5 



solver it is particularly well suited for following discontinu- 
ities in the fluid flow like shocks, contact discontinuities, or 
boundaries between layers of different chemical composition. 
A notable advantage in the present context is its capability of 
solving multi-dimensional problems with high computational 
efliciency and numerical accuracy. Our code makes use of the 
"Consistent Multifluid Advection" (CMA) metiiod (Plewa & 
Miiller 1999) for ensuring an accurate advection of different 
chemical components of the fluid, and switches from the orig- 
inal PPM method to flie HLLE solver of Einfeldt (1988) in the 
vicinity of strong shocks to avoid spurious oscillations (the so- 
called "odd-even decoupling", or "carbunkel", phenomenon) 
when such shocks are aligned with one of the coordinate direc- 
tions in multidimensional simulations (Quirk 1994; Liou 2000; 
Kifonidis et al. 2003; Sutherland et al. 2003). 

Although our hydrodynamic scheme is Newtonian, we have 
included effects of general relativistic (GR) gravity approxi- 
mately in the following way: The gravitational potential used 
in our simulations can be symbolically written as <l)(r, j?) = 
<D^^*t(r, + (Of^(r) - 0^^«"(r)). We compute the Newtonian 
gravitational potential Ojq™' for the two-dimensional axisym- 
metric mass distribution by expanding the integral solution 
of the Poisson equation into a Legendre series, truncated at 
I - 10 (cf. Muller & Steinmetz 1995). General relativistic 
effects are approximately taken into account by the spheri- 
cally symmetric correction term 3)^^ - O^^, where fl)™ de- 
notes an effective general relativistic gravitational potential as 
employed for spherically symmetric simulations (see Rampp 
& Janka 2002, Eq. 53) and O^^"^ is its Newtonian counter- 
part. The general relativistic potential is deduced from 
a comparison of the Newtonian and relativistic equations of 
motion in spherical symmetry and includes terms due to the 
pressure and energy of the stellar medium and neutrinos (see 
Rampp & Janka 2002). Both, (Dfj^(r) and Of™'(r) are com- 
puted using angular averages of the evolved variables. In two- 
dimensional simulations which cover only a limited range of 
latitudes < j^min < ^mia < ^ around the equatorial plane 
we set a>(r,i?) = a>fo(f). 

The source terms Qmt^ 2m,?, 6e> and 2n on the right- 
hand sides of Eqs. (2,3,5,6) are determined by the solution of 
the neutrino transport equations. The source terms Rk depend 
on changes of the composition according to nuclear burning. 
Unless stated otherwise, the EoS we apply is the same as de- 
scribed in detail in (Rampp & Janka 2002, Appendix B). 

Note that PROMETHEUS only solves the left-hand sides 
of the hydrodynamic Eqs. (1-7) and that the EoS is not eval- 
uated during this procedure. The computation of the terms on 
the right-hand sides, i.e the gravitational, neutrino, and burning 
effects, as well as the evaluation of the EoS and, if necessary, 
the determination of the NSE composition, are done in operator 
split steps. 

2.2. Neutrino transport 

Solving the full two-dimensional neutrino transport equation 
would be the most precise way of simulating supernovae in two 
dimensions. On the level of the radiation moments equations 



this would, compared to the ID case, lead to several new de- 
grees of freedom and one additional moments equation, which 
would require five more closure relations to be introduced. 
Therefore we limit ourselves to a less ambitious extension of 
the one-dimensional moments equations to two dimensions in 
which new degrees of freedom (including, in case of the con- 
sidered azimuthal symmetry, the lateral flux and off-diagonal 
pressure tensor terms which account for neutrino viscosity) are 
set to zero, but the additional moments equation (in our case 
for the lateral flux) is taken into account (in the sense that we 
retain the terms corresponding to the lateral neutrino pressure 
gradients at optically thick conditions). In the next section, we 
will justify why we believe this approximation is sufficient for 
simulating two-dimensional neutrino transport in the consid- 
ered models of core-collapse supernovae, but also explain why 
further simpUfications are not allowed. 

The great advantage of our "ray-by-ray" 2D neutrino trans- 
port is that the neutrino moments equations at different latitudes 
(except for some terms which can be accounted for explicitly 
in an operator spfit) decouple from each other. Therefore, for 
each "radial ray", i.e. for all zones of same polar angle, the mo- 
ments equations can be solved independently. Except for some 
additional terms this problem is identical to solving A'^ times 
the moments equations for a spherically symmetric star with 

being the number of grid zones in polar direction. 

We further make the following standard assumptions: First, 
we ignore neutrino oscillations in the SN core, where neutrino 
effects are treated by solving the transport problem. This is jus- 
tified if one ignores the results of LSND; then the parameters 
for atmospheric and solar neutrino oscillations predict the reso- 
nant MSW effect to be at densities far below 10^ g cm so that 
this effect has no influence on the region of interest (for the pos- 
sible importance of neutrino flavor conversions at high densi- 
ties, see Fuller & Qian 2005). Furthermore, non-resonant oscil- 
lations are strongly suppressed in the proto neutron star (PNS) 
due to first and second order refractive effects, see Hannestad 
et al. (2000). Second, the medium, even in the PNS, initially 
does not contain any muons and temperatures and densities 
are always too low to produce tauons, which implies small or 
vanishing chemical potentials for the fi and t type neutrinos. 
Further, the opacities are nearly equal for v^, v^, Vj, and Vt- 
Therefore, we treat these four neutrino types identically and 
set fiy^ = yUy^ = 0. We will notate them collectively as "vx". 



2.2.1. Moments equations 

In a two-dimensional transport scheme assuming azimuthal 
symmetry, the specific intensity lit, r, tp, e, n) does not de- 
pend on the azimuth ip. We describe the direction of propaga- 
tion n by the angle cosine n = n- r/|r|, measured with respect 
to the radius vector r, and the angle cj. Then azimuthal sym- 
metry implies I(. . .,fi,a)) = I(. . .,fi, -ai), see also Appendix 
B. Making the additional assumption that I is independent of 
0), each of the angular moments of the specific intensity can be 
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expressed by one scalar, namely 



ing angular zones. For this purpose we solve the equation 



+1 

{J,H,K,L,...]{t,r,&,e) = ^ J d/i/i'^''^'^ -'JCf, r, e.^*) , 

(8) 

where we have used Eqs. (B.9-B.12) - which follow from our 
assumptions - to reduce the number of independent variables 
in the angular moments of the neutrino intensity as defined in 
Eqs. (B.3). As usual, e denotes the energy of the neutrinos. 
As a consequence of the afore mentioned assumptions the set 
of moments equations for describing the evolution of neutrino 
energy and flux in the comoving frame, given by Eqs. (B.4- 
B.7) in the Newtonian, Oiv/c) approximation, simplifies to 
Eqs. (B.13) and (B.14) for the remaining independent variables 
J and H = Hr. In our approximation Eqs. (B.6) and (B.7) can 
be ignored because the variables the evolution of which they 
describe, i.e. H§ and H^, are strictly set to zero. With !J - Jje, 

= Hje, 'K - Kje, and X - Lie, the moments equa- 
tions describing the evolution of neutrino number are given by 
Eqs. (B.15)and (B.16). 

The velocity-dependent terms, in the order of their ap- 
pearance and grouping in Eqs. (B.13-B.16), correspond to the 
physical efl'ects of radiation advection with the moving stellar 
fluid (expressed by terms containing fird/dr and ij}§lr)dld& 
for the velocity components Vr and respectively), radiation 
compression, Doppler shifting (within the energy derivatives), 
and finally the fluid acceleration. Note that in the current ver- 
sion of our computer code the jSr5/5f-derivatives of the mo- 
ments, which appear at the end of the respective first lines of 
Eqs. (B.13-B.16), cf. Rampp & Janka (2002), are ignored. 

The system of moments equations (B . 1 3-B . 1 6) is very sim- 
ilar to the Newtonian, 0(v/c) moments equations in spherical 
synunetry (see Rampp & Janka 2002, Eqs. 7,8,30,3 1). We have 
set the additional terms arising from our approximative gener- 
alization to two dimensions in boldface. Adding general rela- 
tivistic (GR) efi'ects in the spirit of Rampp & Janka (2002) to 
our approximative 2D transport proves not to alter the boldface 
terms, so that the corresponding equations are Eqs. (54-57) in 
the latter paper, which handle the GR effects for ID transport, 
plus those terms in Eqs. (B.13-B.16) that are typeset in bold- 
face. The equations are closed by substituting K = fx-J and 
L = /l- J, where /k and ft are the variable Eddington factors. 

In order to discretize Eqs. (B.13-B.16), the computational 
domain [0, rmax] X [&mm, ^max] is covered by Nr radial and A^^ 
angular zones, where ■& = Q,n correspond to the polar axis and 
■& = n/2 to the equatorial plane of the spherical grid. All neu- 
trino variables are defined on the angular centres of the zones 
with the coordinate = j(^k + ^^k+i) being defined as 
the arithmetic mean of the corresponding interface values. The 
equations are solved in two operator- split steps corresponding 
to a lateral and a radial sweep. 

In a first step, we treat the boldface terms in the respec- 
tively first lines of Eqs. (B.13-B.16), which describe the lateral 
advection of the neutrinos with the stellar fluid, and thus couple 
the angular moments of the neutrino distribution of neighbour- 



IdE^ 1 5(sinj?j6^5) _^ 
c dt rsind^ && 



(9) 



where H represents one of the moments J, H, J', or 'H. After 
integration over the volume of a zone {i + + |), where i 
and k are the indices of the radial and lateral zones, the finite- 
differenced version of Eq. (9) reads: 



— n+l 

i ,k+ \ 



i+ i ,k+ i 



Ct' 



n+1 



■Ct" 
1 



Ay, 



2 («:+!) 



- aa,,,^,,,,h;,^^^J = o, (10) 



with the volume element Ay,+ i i = y {fi^^ - rf) {cos'&k - 
cos i?t+i) and the surface element AA,-^ i ^ = tt (rf^^ - rf) sin i?^. 
Note that additional indices of H which label energy bins 
and the different types of neutrinos are suppressed for clar- 
ity. To guarantee monotonicity, upwind values of the moments 

l+^,Ki+l/2(k) 



\k-- 



for ^(>,„i^>0, 
else. 



(11) 



are used for computing the lateral fluxes across the interfaces of 
the angular zones. The time-step limit enforced by the Courant- 
Friedrichs-Lewy (CFL) condition is 



A?cFL = min 

i,k 



Ax,,: 



\^«;i+\,k+\ 



(12) 



where Ax,, is the zone width in lateral direction. It turns out that 
this condition is not restrictive in our simulations. In practice, 
the numerical time-step is always limited by other constraints. 

Now, in the second step, the radial sweep for solving 
Eqs. (B.13-B.16) is performed. Considering a radial ray with 
given radial discretization of the equations (and of 

their general relativistic counterparts) proceeds exactly as de- 
tailed in Rampp & Janka (2002). The terms in boldface not yet 
taken into account in the lateral sweep do not couple the the 
neutrino distribution function of neighbouring angular zones 
and thus can be included into the discretization scheme of the 
radial sweep in a straightforward way. 

2.2.2. Eddington factors 

In analogy to the treatment of the moments equations in the 
previous subsection the calculation of the variable Eddington 
factors could in principle also be done on (almost) decoupled 
radial rays. However, such a procedure would require a sizeable 
amount of computer time. On the other hand, the Eddington 
factors are normalized moments of the neutrino phase space 
distribution function and thus, in the absence of persistent 
global deformation of the star, should not show significant vari- 
ation with the angular coordinate (cf. Rampp 2000; Rampp & 
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Janka 2002). Therefore we have decided to determine the vari- 
able Eddington factors only once for an "angularly averaged" 
radial ray instead of computing them for each radial ray sepa- 
rately. The corresponding reduction of the computational load 
can be up to a factor of 10 (Rampp & Janka 2002). 

In analogy to the spherically symmetric case the code 
solves a time-step of the one-dimensional neutrino transport 
on a spherically symmetric image of the stellar background. 
The latter is defined as the angular averages of the struc- 
ture variables ^ £ {p,T, Ye,Pr, • • • } according to ^{t, r) = 
d cos & ^(t, r, ff) wit, r, if), where a - cosi&max) and b = 
cos(j?min) and w is a weighting function which can be either 
p or 1 depending on which choice is more appropriate. The 
computation of the ID transport step proceeds exactly as de- 
scribed for the spherically symmetric version of the code (see 
Rampp & Janka 2002), i.e. the coupled set of the ID moments 
equations and a ID model Boltzmann equation are iterated to 
convergence to obtain solutions for the variable Eddington fac- 
tors /k and /l. These variable Eddington factors are used for all 
latitudes § of the multidimensional transport grid when solving 
the two-dimensional moments equations. We will discuss and 
try to estimate the possible errors associated with our approxi- 
mate 2D transport treatment in Sect. 2.3. 

2.2.3. Coupling to the hydrodynamics 

The system of neutrino transport equations (B. 13-B. 16) is cou- 
pled with the equations of hydrodynamics (1-6) by virtue of the 
source terms 

= y-47rmbysgn(v) de(£f(e), (13) 

y Jo 

/-*CO 

Qe = deCf\e), (14) 

y Jo 

Qur = deC(;,\e), (15) 

y ^ Jo 

A-TT C°° 

Qm» = Tj-— X '^^'^*'v(f)' (16) 

where mby denotes the baryon mass, = e~^Cf\e), and 

Cf\t, r, 1?, e) = (4;r)-i / dQ Cyit, r, §, e, ri) and C^^\t, r, e) = 

(C*.'y',C*/^) = (47r)-' J dQn Cv(f,r,)?,e,K) are angular mo- 
ments of the collision integral of the Boltzmann equation, Cy. 
Note that in Eqs. (13-16) the moments of the colUsion integral 
are summed over all neutrino types v e {ve, Vg, v^,, v^,, v^, v^), 
and sgn(v) = -Hi for neutrinos and - 1 for antineutrinos. 
Remember that we treat v^,, v^,, Vt, and Vr identically because 
their matter interactions are nearly equal. They do not transport 
electron lepton number and therefore do not contribute to 2n- 
In the following we suppress the index v. 

Our simplification of the neutrino transport equations en- 
forces a radial flux vector, i.e. the angular flux component 
Hfi = 0. However, as we shall demonstrate in Sect. 2.3.2, 
the corresponding lateral component of the momentum transfer 
from neutrinos to the stellar medium, described by the source 
term Qm», can not be neglected in the Euler equation of the 



stellar fluid, Eq. (3), when the neutrinos are tightly coupled to 
the medium. This impUes that we should solve the moments 
equation for the lateral transport of neutrino momentum (B.6) 
which, using the assumptions Eqs. (B.9-B.12), simpUfies to 



j8^\„^ 1 d{J-K) 1 dp. 



P» d{J -K) _ d 

2? dt Te^^ 



2\dr r d§ r ' 



(17) 



On the other hand, at the conditions present in the optically 
thick PNS an isotropic neutrino distribution (i.e. J = 3K) and 
diffusion are good approximations. Concerning the velocity de- 
pendent first term in Eq. (17), which is of order ^ j where A is 
the neutrino mean free path, see the discussion about first or- 
der and 2nd order diffusion in Mihalas & Mihalas (1984, §97, 
p. 458ff). Omitting this term means that we ignore all effects of 
neutrino viscosity. Thus the above equation simplifies consid- 
erably to (assuming also stationary conditions, i.e. d/dt = 0) 



deC«(e)= r 
Jo 



de 



1 diJ-K) 
2r dd- 



f 

Jo 



de 



1 dJ 



(18) 



where the terms with the energy derivative in Eq. (17) vanish 
when integrating over the neutrino energy. Given J(r, §) as the 
solution of the moments equations (B.13, B.14), Eq. (18) to- 
gether with the definition Eq. (16) allows us to compute an ap- 
proximation for the momentum exchange rate, Qm§, between 
neutrinos and the stellar fluid. 

Finally, Eq. (B.7) would result in a momentum transfer 
from the neutrinos to the medium in azimuthal direction, C*,'', 
in the presence of rotation, 0. Using our usual assump- 
tions Eqs. (B.9-B.12) we obtain for Eq. (B.7) 



C('>(e) = 



dr 
_ d_ 
d'e 
1 



2c dt ' 2c 



dt 



2c dt 
dr r 



(H-L) 



(19) 



However, again assuming diffusion, stationary conditions, and 
ignoring neutrino viscosity, aU terms of this expression vanish. 



The numerical discretization of Eq. (18) reads 



J. 



J:. 



(20) 



Since Eq. (18) is valid only in the limit of an optically thick 



medium we set C 



(1) 



1+ i ,k+ 



1 = if the density p,+ 1 in a zone 



drops below 10^^ g cm"^. The chosen cut-ofF value is obviously 
specific to the core-collapse supernova problem where outside 
of this density the neutrino pressure gradients, in particular in 
the lateral direction, turn out to be negligibly smaU. 



8 



R. Buras et al.: Two-dimensional hydrodynamic core-collapse supernova simulations witii specb-al neutrino ti'ansport 



2.3. Discussion and tests of ttie numerical sctieme 

In the last subsection we described our implementation of an 
approximative neutrino transport scheme for two-dimensional 
configurations using spherical coordinates and azimuthal sym- 
metry. 

Besides adopting the approximations of general relativis- 
tic efi'ects from the spherically symmetric VERTEX code of 
Rampp & Janka (2002) we made two major approximations 
of the transport equations. First, the dependence of the spe- 
cific intensity on the direction of propagation n is replaced 
by a dependence on only one angle cosine ju. Secondly and 
closely related to the first approximation, we use scalar variable 
Eddington factors. These are obtained from the solution of the 
one-dimensional transport equations on a spherically symmet- 
ric image of the stellar background. Note, however, that our 
treatment described here has been extended considerably by 
additional terms (cf. Sects. 2.2.1, 2.2.3) compared to the sim- 
pler ray-by-ray transport scheme suggested in Rampp & Janka 
(2002, Sect. 3.8). 

In the following we point out limitations of our approach 
and try to critically assess their influence on the results obtained 
with our method. 

2.3.1 . Treatment of general relativity 

The radial neutrino transport contains gravitational redshift 
and time dilation, but ignores the distinction between coordi- 
nate radius and proper radius. This simplification is necessary 
for consistently coupling the transport code to our basically 
Newtonian hydrodynamics. Of course, one would ultimately 
have to work in a genuinely multi-dimensional GR frame- 
work which, among other complications (see e.g. Cardall & 
Mezzacappa 2003), entails abandoning the use of the Lindquist 
metric. 

Tests showed that in spherically symmetric simulations our 
approximations seem to work satisfactorily well (Liebendorfer 
et al. 2005; Marek et al. 2005), at least as long as the infall ve- 
locities do not reach more than 10-20% of the speed of light 
in decisive phases of the evolution. Unless very extreme con- 
ditions are considered (e.g. very rapid rotation) gravity in su- 
pernovae is dominated by radial gradients. We therefore expect 
that efl'ects from a fully multi-dimensional GR treatment are 
small and disregarding them is acceptable in view of the other 
approximations made in the multi-dimensional treatment of the 
transport. Some tests of applicability of our 2D approximation 
of relativistic gravity can be found in Marek et al. (2005). 

It should be noted that the efi'ective relativistic gravitational 
potential which we use in multi-dimensional simulations (see 
Sect. 2.1) cannot guarantee strict conservation of momentum 
(for a detailed explanation see Marek et al. 2005). However, 
in practical applications the size of the violation turns out to 
be very small, see Fig. 1. In a 2D model with lateral periodic 
boundary conditions. Model sl5Gio_32.b, the total momentum 
of the matter on the computational grid (which should remain 
zero) reaches moderate values (corresponding to velocities of 
at most lOOkm/s for a mass of IM©). However, due to the pe- 
riodic boundary conditions the center of mass does not move 



significantly away from the center of the grid because the mat- 
ter which flows out of the grid through one lateral boundary 
reappears at the opposite grid boundary. In a model with lat- 
eral reflecting boundary conditions. Model sll.2Gio_128.b (in 
this model the lateral grid reaches from pole to pole; this model 
will be described in detail in Buras et al. 2005), no such shift 
of mass occurs, so that the center of mass moves according 
to the total momentum. However, as can be seen in Fig. 1, 
for this model the total momentum remains much smaller than 
in Model sl5Gio_32.b (corresponding to velocities of at most 
50km/s for IMq), and the time average of the velocity is close 
to zero. Consequently, the distance between the center of the 
grid and the center of mass remains far below 1km. This be- 
haviour is explained by the fact that the direction of the artifi- 
cial "force" creating the violation of momentum will usually be 
from the center of mass towards the center of the grid. In sum- 
mary, the nonconservation of total momentum due to the em- 
ployed treatment of the effective relativistic gravitational po- 
tential remains a very smaU perturbation which has no influ- 
ence on our simulations. 

2.3.2. Lateral coupling and neutrino pressure 

gradients 

Including the terms printed in boldface in Eqs. (B. 13-B. 16) and 
taking into account the acceleration of the stellar fluid in the 
lateral direction by neutrino pressure gradients, we extend our 
non-equilibrium transport scheme beyond a simple ray-by-ray 
description. This extension is consistent with our fundamental 
assumption of azimuthal symmetry of the intensity I (which 
imphes = 0) and makes sure that our transport description 
produces the correct behaviour in important limiting cases. 

First, in optically thick regions where neutrinos are tightly 
coupled to the stellar fluid, neutrinos must be allowed to be car- 
ried along with (laterally) moving fluid elements. This assures 
the conservation of the total lepton number (Fiep = 7^ + Yy) in 
these fluid elements in the absence of neutrino transport relative 
to the medium (DFiep/ff = for the Lagrangian time deriva- 
tive). Below we will demonstrate that if this conservation is vi- 
olated by omitting the terms describing lateral advection, radia- 
tion compression, and Doppler eff'ects in the neutrino moments 
equations, fluctuations in the lepton number are artificially in- 
duced, which grow and trigger macroscopic fluid motions, e.g., 
by buoyancy-driven rise of high-Fe domains. Secondly, when 
neutrinos yield a significant contribution to the pressure (as is 
the case in the dense interior of the hot, nascent neutron star) 
the inclusion of lateral neutrino pressure gradients is again im- 
portant to prevent artificial acceleration of the fluid by gradi- 
ents of the gas pressure and to provide a restoring force which 
damps the motion of the fluid when it transports neutrinos from 
one place to another. 

In contrast, the omission of angular flux components (H§ = 
0) means the disregard of "active propagation" of neutrinos rel- 
ative to the stellar fluid. Although it is not clear that this is a 
good assumption (for example, the effects of neutrino difiiision 
might be underestimated, see below), it is unlikely to lead to 
fundamental inconsistencies, because it is the correct physical 
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Fig. 1. a Component of the linear momentum in polar axis direction 
of all matter on the computational grid. The solid line shows its time 
evolution for Model s 1 5Gio_32.b, in which a lateral wedge around the 
equator and with periodic boundary conditions was used (this model 
is described in detail in Sect. 3.2.1). The dotted line belongs to Model 
si 1.2Gio_128.b, which is a model with a full 180° grid and thus with 
reflecting boundary conditions at the poles (in lateral direction); this 
model will be described in Buras et al. (2005). b Distance between the 
center of the grid and the center of mass for the same two models. 



limit for situations where the opacity is very high. In the same 
spirit also off-diagional terms of the neutrino pressure tensor 
(Pij with / j) can be dropped, implying that effects of neu- 
trino viscosity are ignored. 

We point out here that another inconsistency is imported 
into our treatment of 2D transport. The omission of angular flux 
components causes the problem that the correct limit at large 
radii and small optical depth may not be accurately reproduced 
because the evolution of the radiation moments J and H even 
at large radius depends on lateral gradient terms that include 
fy. Fortunately, the specifics of the supernova problem help us 
justifying our approach: At large radii and low optical depth 
the lateral component of the fluid velocity is usually small 
(in particular relative to the radial component of the velocity, 
which determines relevant timescales) so that all terms scaling 
with - Vfi/c become of minor importance. Secondly, multi- 
dimensional effects of neutrino transport usually lead to a neu- 



trino distribution outside of the neutrinosphere, which is more 
isotropic than the stellar structure (Livne et al. 2004), an ef- 
fect which our transport treatment tends to underestimate any- 
way. Thirdly, at large radii the non-cancelling yS^-terms are sup- 
pressed by a factor 1/r, and the small remaining lateral redis- 
tribution of the neutrino flux is irrelevant because the couphng 
between neutrinos and stellar medium (e.g. neutrino heating 
and cooling) vanishes. 

Criterion for instability In the absense of neutrino diffusion the 
hydrodynamic stability in the neutrino trapping regime of the 
PNS is tested by the Ledoux criterion. 



\ds jyp dr \oY i^p dr 



(21) 



with s - 
tropy Sy 



s + Sy being the entropy including the neutrino en- 
and Y - Fiep - Y^^ + Yy being the total lepton num- 
ber. Stability means Cl < 0. For 2D models Eq. (21) is evalu- 
ated on a spherically symmetric image of the star as defined in 
Sect. 2.2.2. The corresponding Brunt- Vaisala frequency 



a>BV = sign(CL) 



|Cl 



(22) 



where g = -dO/dr is the gravitational acceleration, is closely 
related to Cl and denotes the growth rate of fluctuations, if pos- 
itive (unstable modes), and the negative of the oscillation fre- 
quency, if negative (stable modes). 

Recently, Bruenn et al. (2004) presented a more elaborate 
discussion of hydrodynamic stability in the PNS including the 
effects of neutrino diffusion (an extension of a previous anal- 
ysis by Bruenn & Dineva 1996). They argue that local pertur- 
bations in the lepton number will be reflected in the neutrino 
phase space and thus cause a net neutrino diffusion which tries 
to wash out the perturbation, an effect which can be accounted 
for by a "response function". Since neutrinos also carry entropy 
the neutrino diffiision that smoothes the lepton number pertur- 
bation will create an entropy perturbation. This effect is charac- 
terized by a "cross response function". Of course, entropy per- 
turbations will analogously induce an equiUbrating net neutrino 
diffusion which at the same time carries lepton number, corre- 
sponding to another "response function" and a "cross response 
function". Bruenn et al. (2004) found in a numerical analysis 
that perturbation-induced neutrino diffusion transports lepton 
number more efficiently than entropy, and that the transport 
of lepton number reacts faster to entropy perturbations than to 
lepton number perturbations. For such a situation convective 
instability should set in at most stellar conditions, even when 
the fluid is Ledoux stable. In particular, Bruenn et al. (2004) 
describe two kinds of instabilities in the presence of neutrino 
diffusion, "lepto-entropy finger" (LEF) convection and "lepto- 
entropy semi-convection" (LESC). 

Unfortunately, Bruenn et al. (2004) did not provide detailed 
information about the values of the (cross) response functions 
in the different regions of the PNS, and the calculation of these 
response functions appears to be quite involved. Although one 
might use their approximate values for the response functions 
as an estimate, Bruerm et al. (2004) show that these depend 



10 



R. Buras et al.: Two-dimensional hydrodynamic core-collapse supernova simulations witii specb-al neutiino ti'ansport 



on which interaction rates are taken into account. They further 
mention that the existence of LESC is sensitive to the exact val- 
ues of the cross response functions. Due to these uncertainties 
we refrain from applying their stabiUty analysis here. Instead, 
in order to account for efficient transport of lepton number hap- 
pening in convective 2D models between buoyant fluid ele- 
ments and their surroundings, we extend Eq. (21) to a "Quasi- 
Ledoux" criterion, following Wilson & Mayle (1993) and Keil 
(1997): 



Cql 



-At) 



d{s} 



dp 

dY 



d{Y} 
dr 



Here, the brackets {) denote angular averaging to emphasize 
that the last term, dV/dr, is evaluated locally. The additional 
parameter j8 has been introduced as a measure of the efiiciency 
of the exchange of lepton number between a perturbed fluid 
element and its surroundings. From simulations it was found 
that values of yS ^ 0.5-1.5 are typical, and that the criterion is 
not very sensitive to the exact value of j8 (Keil 1997). Therefore, 
we will use yS = 1. 

Difl'erent from convection inside the PNS (i.e. below the 
neutrinosphere) exchange of entropy and lepton number by 
neutrino transport between moving fluid elements and their sur- 
roundings does not play an important role in the gain layer, 
where neutrinos and stellar fluid are much less tightly coupled. 
Therefore we apply Eq. (21) using the medium entropy s and 
F = Fe when evaluating the stabiUty criterion, or the Brunt- 
Vaisala frequency, at densities lower than lO'^ gcm"^, con- 
sidering the latter as the density below which neutrinos begin 
decoupUng. 



0.05 



0.04 



0.03 



0.02 r 



0.01 



0.00 



a 




Stability analysis In order to assess the importance of the lat- 
eral gradient terms which we take into account in our neutrino- 
hydrodynamics treatment, we performed three test calculations 
which were started by imposing random seed perturbations in 
the density with a maximum amplitude of +2.5% on an early 
post-bounce model which was taken from a spherically sym- 
metric simulation. The perturbed model was then evolved for 
a few milliseconds (corresponding to a multiple of the relevant 
dynamical timescale) in order to test diflBrent variants of the 
2D transport code described above. 

Figure 2a shows the standard deviation of density fluctua- 
tions in lateral direction. 



Pk(r) - {p(r)},.) 



(24) 



where Aci?^. = cosid-f^ \_) - cos(i9'j^i). A' = J^^, Acj?/.-, and 

ipir))^ = Yi^, pk'{r)Ac§k', as a function of radius after 3.6 
ms of 2D evolution. The quantity cTp serves as a convenient 
measure to specify the magnitude of the density fluctuations. 
Employing the full implementation of our neutrino treatment, 
which takes into account both lateral gradients in the neutrino 
moments equations as well as lateral gradients of the neutrino 
pressure in the fluid equations, the initial perturbations do not 
grow to large instabilities anywhere inside the PNS, which is 



Fig. 2. a Standard deviation of the density a-p (see Eq. 24) indicating 
convective activity inside the neutron star. The ID Model sl5Gio_ld.b 
was mapped to 2D (16 zones with a resolution of 2.7°) 27ms after core 
bounce and the density distribution p was perturbed. The plot shows 
the situation after 3.6ms of dynamical evolution computed with differ- 
ent implementations of the 2D transport equations. The dotted line is 
the initial value of cr^, the thick solid line shows the standard deviation 
(Tp when the lateral terms are included in our scheme as described in 
Sects. 2.2.1 and 2.2.3. For comparison, the thin solid line shows CTp 
when running with pure ray-by-ray transport (i.e., without all bold- 
face terms in Eqs. (B.13-B.16) and without the lateral component of 
neutrino pressure gradients), and the dash-dotted Une corresponds to 
a simulation where the lateral terms of Sect. 2.2.1 and Appendix B 
were taken into account but not the neutrino momentum transfer to 
the fluid discussed in Sect. 2.2.3. b Brunt- Vaisala frequency for the 
same model at the beginning of the test calculations, derived from 
the Quasi-Ledoux criterion, Eq. (23), using s = s + Sy and Y = Fiep. 
Negative wbv indicates convectively stable regions. 



in accordance with the prediction by our Quasi-Ledoux crite- 
rion. The residual fluctuations which can be discerned at radii 
between 15km and 27km seem to be of physical origin because 
the Quasi-Ledoux criterion predicts instability in part of this 
region, see Fig. 2b. 

When switching off the effects of the lateral component of 
the neutrino pressure gradient we notice a strong amplifica- 
tion of the initial perturbations in a region where the damp- 
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ing should be strongest according to the Quasi-Ledoux crite- 
rion. A naive ray-by-ray scheme, which in addition to lateral 
neutrino pressure differences also disregards the lateral gradi- 
ents in the neutrino moments equations (including those terms 
which account for the lateral advection of neutrinos), produces 
spurious convective activity in a broad region between 15km 
and 25km. The latter phenomenon is not Ledoux convection 
because it transports lepton number and entropy in the wrong 
direction (i.e. against the gradients). While lateral fluid motions 
(adequately modelled by the two-dimensional hydrodynamics 
scheme) tend to damp fluctuations, the resulting advection of 
electron number produces artificial flucuations of lepton num- 
ber and therefore obviously instigates buoyancy motions, if the 
advection of electrons is not accompanied by the correspond- 
ing advection of neutrinos (i.e. if the total lepton number is 
not preserved in the moving fluid elements even when diffu- 
sion of neutrinos is unimportant). This is the case when the 
terms accounting for lateral advection in the neutrino moments 
equations are neglected. Neutrino pressure differences, on the 
other hand, produce restoring forces which help damping lat- 
eral motions after local fluctuations have been reduced by fluid 
advection. 

We have not attempted to perform a very detailed analysis 
of all efTects of lateral advection and pressure gradients of neu- 
trinos to elaborate our understanding beyond the more qualita- 
tive insights described above. However, we interprete our tests 
as a demonstration that approximations of multi-dimensional 
transport schemes must be tested carefully for the possibility 
of producing spurious convective activity (or suppression of 
convection) in the newborn neutron star, where neutrinos con- 
tribute significantly to the total pressure and total lepton num- 
ber density. 

2.3.3. Lateral propagation of neutrinos 

A number of our two-dimensional supernova simulations, in 
particular those which produced lively hot-bubble convection, 
showed transient neutrino-bursts when narrow downflows of 
accreted stellar gas entered the cooling region and penetrated 
down to the vicinity of the neutrinosphere. In our simulations 
such bursts occurred roughly every 20 ms, typically persisting 
for a few ms. The neutrinos of such bursts, which are the result 
of locally enhanced neutrino emission, would naturally prop- 
agate in all directions in a fully multi-dimensional treatment 
and would therefore illuminate the surroundings in all direc- 
tions. Our code ignores the lateral propagation of these neu- 
trinos, thus the burst is radiated away in the radially outward 
direction and its lateral width is essentially constrained to the 
layers above the hot, radiating area. Here we shall argue that al- 
though locally and transiently the neutrino heating rates can be 
incorrect by up to a factor of two (in very rare, extreme cases) 
due to the disregard of lateral neutrino propagation, the transfer 
of energy between neutrinos and stellar gas is not significantly 
changed on larger spatial and temporal scales, and therefore the 
global dynamics of our supernova simulations is not likely to 
be affected significantly. 



Basically, our simplified neutrino transport overestimates 
the heating in the radial direction at polar latitudes where "hot 
spots" occur, while adjacent "rays" experience less heating 
than in a 2D neutrino transport which includes lateral neu- 
trino fluxes. Truly two-dimensional transport tends to redis- 
tribute neutrinos in lateral direction, in particular in the semi- 
transparent and transparent regimes where the mean free path 
becomes large. This can lead to a more uniform spatial distribu- 
tion of neutrinos exterior to the neutrinosphere than in case of 
our ray-by-ray treatment, in particular in the presence of local 
hot spots (see also the discussion in Livne et al. 2004). 

To test the implications of such a lateral neutrino redistribu- 
tion for the neutrino heating behind the shock, we performed a 
post-processing analysis in which we averaged the frequency- 
dependent neutrino densities over all polar angles at a given 
radius and time and used the result to recalculate the local net 
heating rates. Then we compared the recalculated heating rate 
with the actual heating rate of the simulation, both integrated 
over the respective gain layer, for different times. Hereby the 
gain layer was defined for each grid value of ■& separately as 
the region between the shock and the gain radius, i.e. the inner- 
most radial point where net heating occurs. For this purpose, 
the position of the gain radius was also redetermined for the 
recalculated heating rates. 

We have carried out this analysis for an 11 .2Mq progen- 
itor (si 1.2, Woosley et al. 2002) computed with 32 angular 
zones and ~ 90° lateral wedge during its post-bounce evolu- 
tion (Buras et al. 2003b), a model which showed powerful hot 
bubble convection. The evaluation was started at fpb = 100ms; 
before that time hot bubble convection was weak and no strong 
local bursts of accretion luminosity happened. We found that 
the angular averaging of the neutrino densities hardly changes 
the total net heating rate in the gain layer, 5f£gi, see Fig. 3, up- 
per panel. At very few instants of evaluation we discovered a 
significant increase of 5iE„i by at most 30%. In the temporal 
average this difference is reduced to only a few percent. For the 
average net heating per baryon in the gain layer (Fig. 3, middle 
panel), which, concerning global supernova dynamics, we con- 
sider to be the more decisive quantity (Janka 2001), the values 
for the two heating rates are almost indistinguishable. 

Interestingly, the differences are somewhat larger if we sep- 
arately analyze downflows of cold material and hot bubbles in 
the gain layer; a zone of the numerical grid is attributed to a 
downflow if the negative velocity is more than 1.5 times the 
angle-averaged (negative) velocity (v)^ at a given radius r, oth- 
erwise we define the zone to belong to a high-entropy bubble. 
For the downflows (see also Fig. 3, upper panel, lower curve), 
the time-averaged total net heating decreases by 14% in case 
of the angle-averaged neutrino density. On the other hand it in- 
creases by 12% for the high-entropy bubbles. 

The results from this analysis should be read carefully: It 
does not take into account the dynamical consequences of the 
altered heating rates. Also, the angle -averaging of the neutrino 
density overestimates the spreading of the burst neutrinos in all 
directions, especially close to the gain radius, where we obtain 
by far most of the heating. The nevertheless modest sensitivity 
of the integral quantities to the described averaging of the neu- 
trino distribution, however, gives us some confidence that our 
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Fig. 3.The upper panel compares total (volume integrated) 
neutrino heating rates in the gain layer <5,£gi. The bold up- 
per solid line shows the evolution of 6,E^i for the 11.2Mo 
supernova model of Buras et al. (2003b). The post-processed 
heating rate obtained by averaging the neutrino densities 
over all latitudes is represented by diamond-like symbols. 
For a consistency check crosses are drawn for the post- 
processed neutrino heating rate without using the lateral av- 
erages of the specific neutrino density. Ideally the latter post- 
processed results should be equal to the values returned from 
the model run. Remaining differences result from the fact 
that the post-processing may yield a gain radius shifted by 
one radial zone. For drawing the second curve (thin solid 
line) and the corresponding symbols the analysis was re- 
stricted to downflows in the gain layer, i.e. for regions with 
velocity v < 1.5 (v)^ < 0. The middle panel shows the aver- 
age net heating rate per baryon (the symbols have the same 
meaning as in the upper panel), i.e. the total neutrino heat- 
ing rate in the gain layer divided by the total mass Mgi con- 
tained in this region. The lower panel shows the maximum 
neutrino flux as measured slightly above the gain radius nor- 
malized to the flux average over all latitudes at the same 
radius, = max^ {F) / <F)^. The values indicate the rel- 
ative strength of localized luminosity outbursts. 



approximation of the moments equations is a reasonable step 
towards fully consistent multidimensional simulations. 

Even more difficult to assess quantitatively is the effect of 
replacing the Eddington tensors (c.f. Appendix B) by scalar 
variable Eddington factors which are calculated from a solu- 
tion of the one-dimensional transport equations using a spher- 
ically symmetric image of the star. We have to rely on the 
fact that these variables are normalized moments of the neu- 
trino phase space distribution and thus, in the absence of per- 
sistent global deformation of the star, should not show signif- 
icant variation with the angular coordinate (cf. Rampp 2000; 
Rampp & Janka 2002). Anyway, as we have seen, ignoring 
the lateral propagation of neutrinos does not seem to largely 
affect the neutrino-matter coupling in the semi-transparent re- 
gion. Further, Thompson et al. (2005) have recently shown that 
neutrino viscocity, which is connected with the non-diagonal 
elements of the Eddington tensors, is of minor importance. 
Therefore we hope that the transport is also sufficiently insen- 
sitive to our approximation concerning the variable Eddington 
factors. Final answers can only be expected from simulations 
employing a fuUy multidimensional treatment of the neutrino 
transport. 

2.4. Equation of state 

All calculations presented in this paper were run with the equa- 
tion of state (EoS) described in detail in Rampp & Janka (2002, 
Appendix B). The EoS of Lattimer & Swesty (1991) (LS EoS) 
is employed to treat matter at densities above 6 x 10^ gem"-'. 
Recently it became clear that this EoS underestimates the frac- 
tion of alpha particles (due to an error in the definition of the 
alpha particle rest mass; J. Lattimer, personal communication) 
and concerns were expressed that deficiencies of this EoS in the 
density regime between 10^ g cm"^ and 10^^ g cm~^ could have 



consequences for the possibility of an explosion, which we fail 
to get with the LS EoS (C. Fryer, personal communication). 

In order to test this we have implemented a new EoS in our 
code. This EoS is similar to the low-density EoS (Janka 1999) 
described by Rampp & Janka (2002, Appendix B), but the very 
approximative handling of nuclear dissociation and recombina- 
tion in nuclear statistical equilibrium (NSE) in that EoS is now 
replaced by a composition table in the (p, T, 7e)-space. This ta- 
ble was computed assuming matter to be in NSE and solving 
the Saha equations for equilibrium of a mixture of free neutrons 
and protons, "^He, and ^'*Mn as a representative heavy nucleus 
(Janka 1991). We use the table for T > 5 x lO'K, where NSE is 
a reasonably good assumption. Given the nuclear composition, 
variables such as pressure and entropy can be calculated from 
the EoS of Janka (1999). 

A comparison of this simple 4-species NSE table with a 
sophisticated NSE solver that takes into account 32 different 
species of heavy nuclei showed excellent agreement in entropy, 
pressure, and helium mass fraction for typical conditions met 
in the post-shock layer. Note that the 4-species table overesti- 
mates the mass fraction of free nucleons in parameter regions 
where heavy nuclei dominate the composition and Fg is far 
from 25/54 ^ 0.463 (e.g. for T = IMeV, = 0.3 this is the 
case atp>3xl0'"g cm \ see Fig. 4). This problem originates 
from the restriction to one fixed, representative heavy nucleus 
and the necessity to fulfill charge neutrality. However, NSE in 
the post-shock layer of a supernova is characterized by condi- 
tions where essentially all heavy nuclei are dissociated into free 
neutrons, protons, and a-particles. 

Comparing our composition table with data from the LS 
EoS (Fig. 4) we indeed observe a significantly larger helium 
mass fraction for certain combinations of density and temper- 
ature. Moreover, the region where helium contributes signifi- 
cantly to the composition is larger. This also affects the pres- 
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Fig. 4. Comparison of the LS EoS (thick lines) with the four-species 
NSE table introduced here (thin lines), for = 0.3 and T = 1 (solid), 
2 (dotted), and 3 (dashed) MeV. X„ and Xj, are the mass fractions of 
a-particles and of the representative heavy nucleus, respectively. 



sure and entropy, especially at low temperatures and entropies 
(Fig. 4). 

To test whether the sizeable differences between the EoS 
using our NSE table at densities below 10^' g cm and the LS 
EoS have any influence on supernova simulations we ran the 
spherically symmetric 11.2M0 model of Buras et al. (2003b) 
with both equations of state. This comparison directly tests the 
influence of the different EoS implementations on the dynami- 
cal evolution of the supernova. 

Surprisingly, both runs showed virtually no difference in 
the post-bounce history of the supernova (Fig. 5) and only 
small differences in the post-shock structure (Figs. 6 and 7), 
even though a post-processing analysis where values of T , p, 
Fe were crosswise fed into the two variants of the EoS revealed 




20 40 



60 80 
tpb [ms] 
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Fig. 5. Shock positions for the two dynamical simulations of the post- 
bounce evolution of an 11.2Mo star with the LS EoS (dashed) and 
the four-species NSE table (solid) used in the density regime below 
10" gcm-^ 



differences in entropy and pressure of up to 15% and 10%, re- 
spectively, below the shock (Fig. 8). This phenomenon has a 
simple explanation: On the one hand the relatively high en- 
tropies behind the shock in our models imply that the influence 
of or-particles is in general moderate (Fig. 4). On the other hand 
the accretion layer behind the supernova shock can be consid- 
ered as an approximately hydrostatic situation where the veloc- 
ities of the shocked medium and of the shock itself are much 
smaller than the sound speed c^. In this case the pressure profile 
behind the shock is determined by constraining boundary con- 
ditions, i.e. the mass and radius of the proto neutron star (and 
therefore its gravitational field) on the one hand, and the mass 
accretion rate at the shock on the other (Janka 2001). Keeping 
in mind that the Rankine-Hugoniot conditions coimect the flow 
in the infall region with the conditions in the post-shock layer, 
the pressure profile in the accretion layer is fully determined. 

Thus, for a given pressure profile other variables, such as p 
and T, may vary for different EoSs (see Fig. 7) without chang- 
ing the dynamical evolution of the supernova. Only if these 
variables have an immediate influence on the evolution, e.g. by 
temperature-dependent neutrino emission, do we have to worry 
about the differences and consequences of the composition. 
Neutrino heating, however, is dynamically not very important 
in the gain layer of our non-exploding models, while in the 
cooling region the composition differences between the EoSs 
are negligibly small. 

We conclude that, although the LS EoS has problems with 
the a-mass fraction and other dependent variables in the low- 
density regime at supernova conditions, we do not find a major 
effect on the dynamics of our simulations. 

o 



2.5. Entropy wiggles 

All the models presented here feature a numerical problem, the 
origin of which we found only recently. In all profiles of ID 
calculations, and also partly in 2D calculations (when convec- 
tive effects did not wash out the effect), artificial wiggles ap- 
pear behind the shock. These wiggles are most pronounced in 
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Fig. 6. Comparison of the profiles for the two dynamical simulations 
of the post-bounce evolution of an 1 1 .2Mq star with the LS EoS (dot- 
ted) and the four-species NSE table (solid) used in the density regime 
below 10" gcm"^ at two post-bounce times. 



the entropy and pressure profiles with amplitudes of up to 10%, 
see e.g. Fig. 8. 

The problem is connected with the use of the specific en- 
ergy s - e + ekin in the hydrodynamics equations (1-7), which 
are solved in the PROMETHEUS code with a Riemann solver 
technique, see Sect.^.l, and the corresponding calculation of 
the required index Fg = ^ + 1 with an arbitrarily normal- 
ized energy e instead of the internal energy eint. The Riemann 
solver expects F^ to be calculated with the internal energy Cmt, 




r [km] 

Fig. 7. Relative difference of temperature and density between the sim- 
ulations described in Fig. 6 at ? = 50ms post-bounce. Ar > means 
higher T in the simulation with the four-species NSE table (soUd line 
in Fig. 6) used in the density regime below 10" g cm"^, dito for Ap. 
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Fig. 8. Enlarged comparison of the entropy and the pressure (renor- 

malized with r*) between the two simulations shown in Fig. 6 at 
t = 50ms. The thin lines correspond to the actual values, the thick 
lines correspond to a post-processing analysis of the profiles using the 
respectively other variant of the EoS. 



because for thermodynamical consistency between F^ and the 
adiabatic index Fad = (d In p/d Inp)^, the energy e and the pres- 
sure p should go to zero simultaneously. On the other hand, the 
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Fig. 9. a Shock trajectories of Model sl5Gio_ld.b, described in 
Sect. 3.1.1 (solid), and a calculation with identical micro-physics but 
with our new definition of the energy e (dashed). Note that the visible 
very small differences actually originate from slight differences of the 
gravitational potentials used in both simulations, b Entropy profiles 
for Model sl5Gio_ld.b (thin) and the simulation with our new defini- 
tion of the energy e (thick) at the times 20ms, 40ms, and 60ms after 
bounce (from bottom to top). 

physics should only depend on energy differences, and the en- 
ergy normalization should not play a role, i.e., the EoS as well 
as the hydrodynamics should remain unchanged, independent 
of the chosen zero point of the energy. Therefore it is allowed 
to add an arbitrary normalization to e, still retaining the validity 
of Eqs. (1-7). In our code, we use the definition 

e - eint + erm + CO , (25) 

where erm is the specific baryon and electron rest mass energy 
and eo is a constant offset chosen to be in accordance with the 
energy definition in our high-density EoS provided by Lattimer 
&Swesty (1991). 

r,, as defined above is therefore different from the thermo- 
dynamical index r^, = + 1, which is expected to be fed into 
the Riemann solver. In particular, in regions with relatively low 
internal energy and large value of (erm + eo) (e.g., just behind 
the shock where the matter consists mostly of free nucleons), 
te with our chosen normalization of e can have values which 
are significantly below the physically meaningful lower bound 
of 4/3. Empirically, we found that the numerical noise appears 
whenever < 4/3, which turned out to be the case in the 
post-shock region. However, we did not make efforts to locate 
the source of these numerical fluctuations more precisely than 
linking it to the use of a Riemann solver technique in cases 
where adopts unphysical values. 
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Fig. 10. Progenitor structure. Temperature T, entropy per baryon .v, 
electron fraction Y^, density p, radius r, and shell binding energy iJ^M > 
defined in Eq. (26), as functions of the enclosed mass. At 1.42 M© the 
interface between the silicon shell and the oxygen-rich silicon shell 
is located. T, Y^, and p are taken from model sl5s7b2 of Woosley & 
Weaver (1995), while the s and are computed by applying our 
EoS. 



We have solved this problem by renormalizing the energy 
used in the integration scheme of the hydrodynamics equations 
relative to e by subtracting from e the additional terms e^m + co 
of Eq. (25). The numerical procedure is described in detail in 
Appendix C. In a test calculation with this improved handling, 
the wiggles turned out to disappear, see Fig. 9b. 

Figure 9 also shows that the dynamics and global evolution 
of the models are not affected by the wiggles. This confirms 
that the hydrodynamics solutions given by PROMETHEUS are 
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not globally flawed even if e from Eq. (25) and deviate sig- 
nificantly from each other. Therefore the simulations presented 
in this paper are sound despite the wiggles in the postshock 
profiles (see also Liebendorfer et al. 2005 for a confirmation). 
Of course, future calculations will be performed using the im- 
proved treatment of Appendix C. 



3. Results 

Here we present results from one-dimensional and from first 
two-dimensional simulations with our "coupled ray-by-ray" 
approximation of spectral neutrino transport and with an im- 
plementation of neutrino opacities which is improved com- 
pared to Bruenn (1985). First, we will discuss the differences 
resulting from our enlarged set of opacities in the context of 
spherically symmetric (ID) models. We will also discuss the 
influence of general relativistic corrections on the core col- 
lapse and post-bounce evolution and will investigate the im- 
portance of velocity-dependent terms in the neutrino momen- 
tum equation, especially for a quantitatively correct descrip- 
tion of neutrino heating and cooling between proto neutron star 
(PNS) and supernova shock. In this context we shall present 
an interesting ID model which shows large-amplitude oscil- 
lations of the shock position and neutrino luminosities. Then 
we will elaborate on our first two-dimensional (2D) neutrino- 
hydrodynamical simulation with spectral neutrino transport. 
Finally, we present a 2D model that develops an explosion 
when we retreat from the most complete version of our neu- 
trino transport treatment by neglecting the velocity-dependent 
terms in the neutrino momentum equation. 

Progenitor and model notation: The names of our mod- 
els are chosen with the aim to provide information about 
the employed input physics. All models in this paper were 
started from the progenitor "sl5s7b2", a star with a main- 
sequence mass of 15 Mq kindly provided to us by S.Woosley 
(Woosley & Weaver 1995), see Fig. 10. Thus our model names 
start with "sl5". The simulations were performed either with 
Newtonian gravity ("N") or with our approximative implemen- 
tation of general relativistic gravity ("G"). While most simu- 
lations included the most advanced description of neutrino in- 
teractions ("io"), we performed a number of calculations with 
the standard opacities of Bruenn (1985) plus nucleon-nucleon 
bremsstrahlung ("so"). One-dimensional models are labelled 
with "_ld", the names of the two-dimensional models give in- 
formation about the number of lateral zones N-a (in this paper 
only "_32"). Combining this with the chosen size of the angu- 
lar wedge allows one to infer the equidistant angular zoning of 
the model. Finally, it is very important whether a model em- 
ployed our full implementation of the neutrino transport. Case 
B (".b"), or whether the radial velocity-dependent terms (ex- 
cept for the advection terms) in the neutrino momentum equa- 
tion were omitted. Case A (".a"). The models presented in this 
paper are listed in Table 1. The collapse phase was calculated 
with a Lagrangian grid, the post-bounce phase with a Eulerian 
grid, and the simulations were run with 400 radial hydro zones 
and 230 radial neutrino transport zones, both with geometric 



spacing. The number of tangent rays employed was 250. Inside 
of 400km, the hydro and transport zones were chosen to be 
identical, outside of 400km only a few transport zones were 
added to correctly simulate the neutrino propagation out to the 
outer boundary of the grid at lO'^km. All models were run with 
17 geometrically spaced neutrino energy bins between 2MeV 
and 380MeV. When necessary, the radial grids were refined to 
account for the steepening density gradient at the surface of the 
PNS. 

Definitions: We define the shock radius rghock as the radial 
position where the radial velocity in the shock "discontinu- 
ity" is half of the minimum preshock value of the velocity. 
The gain radius r^^m is the zone interface below the innermost 
radial zone (at given value of v) outside of which no zones 
with net cooling exist. Thus the gain layer is between rgain and 
''shock while the cooling region is below r^^m and includes the 
PNS. The total net energy transfer rates between neutrinos and 
medium, ^fiicooi and 5,iigi in the cooling region and gain layer 
are defined as integrals over the respective regions. An interest- 
ing quantity is the "shell binding energy" E^^^, i-e. the energy 
needed to lift all material above a considered radius out of the 
gravitational potential of the mass enclosed by this radius and 
thus move the material to infinity. In the approximation with 
Newtonian ID gravity (which is a very good assumption for 
the progenitor) its definition is 



gshell/ = 



An 



COS 1?m 

ax - cos 



/-•OO --COS ■ 

Jr Jcosfi, 



4m ('■'' ^Mr', ^ COS 1? d/ , (26) 



where ej*^^' is our so-called "local specific binding energy", 

= ^int(n'?) + 5(v? + v^+v2)(r,^) + (I)-^'°-V). (27) 

The gravitational potential <l)j™'°^'^''(r) is calculated taking into 
account only the mass inside of the radius r and assuming it to 
be distributed spherically symmetrical. These energies, as well 
as the neutrino luminosities, are generally given in the unit of 
FOE, the abbreviation for "ten to the fifty-one erg". Finally, the 
spectrally averaged transport optical depth is evaluated accord- 
ing to the formula 

A;\r', e)H,{r', e)de / H,{r' , e)de , (28) 

where X~^{r',e) = Y^iK] inverse mean free path for 

neutrinos of energy e and type v including all scattering and 
absorption processes. The different contributions to the sum 
are taken from the RHS of the neutrino momentum equation 
as implemented in our transport code. We define an "energy- 
averaged" neutrinosphere (radius) by the condition Tv(ry) = 1. 
Where needed, the energy-dependent neutrinosphere is defined 
by Ty(rv(e), e) = /l^'(r', e)dr' = 1 . Many properties of 2D 
models are evaluated from ID radial profiles, obtained by lat- 
erally averaging the 2D data, see Sect. 2.2.2 for the definition 
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Table 1. Input physics for the sample of computed models presented here. 



Model 


Dim. 


Gravity 


V Reactions 


Transport 


Wedge" 


&-zones 


sl5Nso.ld.b 


ID 


Newtonian 


standard 


CaseB 






sl5Gso_ld.b 


ID 


approx. GR 


standard 


CaseB 






sl5Gso_ld.b* 


ID 


approx. GR 


standard'' 


CaseB 






sl5Gio_ld.a 


ID 


approx. GR 


improved'^ 


Case A 






sl5Gio.32.a 


2D 


approx. GR 


improved 


Case A 


±43.2egr 


32 


sl5Gio_ld.b 


ID 


approx. GR 


improved 


CaseB 






sl5Gio_32.b 


2D 


approx. GR 


improved 


CaseB 


+43.2egr 


32 



^ Angular wedge of the spherical coordinate grid around the equatorial plane. 
Calculation without neutrino-pair creation by nucleon-nucleon bremsstrahlung. 
Calculation without the neutrino-antineutrino processes of Buras et al. (2003a). 
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Fig. Il.Trajectories of mass shells for Model 
sl5Gio_ld.b. The plot also shows the neutri- 
nospheres for (thick solid line), Ve (thick dashed), 
and Vx (thick dash-dotted), the mass shell at which 
the silicon shell becomes oxygen-rich (knotted solid 
line, at 1 A2Mq), and the shock (thick solid line with 
superimposed dashes). Further we have marked the 
regions with a mass fraction of more than 60% in 
iron-group elements (dark shaded). Also shown are 
regions with a mass fraction between 30% and 60% 
in ''He (light shaded). Note that the white region at 
fpb > 200ms and r ~ 90km corresponds to < 0.6 
and < 0.3. Finally, the lower thin dashed Une 
marks the gain radius, while the upper one marks the 
interface between our high-density and low-density 
EoSs (i.e. the upper thin dashed line corresponds to 
a density of p = 6 x 10' g cm"^). 



of this procedure. This also includes the definition of "mass 
shells" in 2D, i.e. a mass shell corresponds to the radius of the 
sphere which encloses this mass. For 2D simulations, the crite- 
rion of convective instability (see Sect. 2.3.2) is also evaluated 
with laterally averaged variables. 

3. 1. Spherically symmetric models 
3.1 .1 . The reference model 

This one-dimensional model, sl5Gio_ld.b, was run with our 
most elaborate "state-of-the-art" description of neutrino opaci- 
ties described in Appendix A and references therein. This treat- 
ment of neutrino-matter interactions is improved compared to 
the rates of Bruenn (1985), which we consider as the "stan- 
dard case", by including effects of nucleon recoil and thermal 
motions, weak magnetism corrections (Horowitz 2002), and 
nucleon-nucleon correlations (Burrows & Sawyer 1998, 1999) 
in neutrino-nucleon interactions. Moreover, the reduction of the 
effective nucleon mass at nuclear densities (Reddy et al. 1999) 
and the quenching of the axial-vector coupling in nuclear mat- 
ter (Carter & Prakash 2002) are taken into account. In addi- 
tion, we include nucleon-nucleon bremsstrahlung and neutrino- 



neutrino interactions Uke scattering as well as pair creation pro- 
cesses between neutrinos of different flavors. The full list of 
included rates and corresponding references can be found in 
Table A. 1. 

We applied the approximations of general relativity de- 
scribed in Sects. 2.1 and 2.3.1 (see also Rampp & Janka 2002 
and Marek et al. 2005). Finally, the simulation was done includ- 
ing the velocity-dependent terms in the neutrino momentum 
equation, i.e. the terms proportional to /3r in Eqs. (B.14,B.16). 
Despite of being formally small (fir < 0.1) these terms turn out 
to be important in supernova simulations as we wiU elaborate 
in Sect. 3.1.3. 

The evolution of this model can be seen in Fig. 1 1 . The 
progenitor needs 174ms to reach core bounce. The typical time 
evolution of the central values of Fiep and s during collapse can 
be seen as functions of the central density in Fig. 12. The shock 
is created at an enclosed mass of M^c = 0A9Mq and a radius 
of = 10.6km. Its prompt expansion stalls after about 1ms at 
r = 32km, turning the shock into an accretion front. Initially 
the high mass accretion rate through the shock leads to matter 
being accumulated between the PNS and the shock; the energy 
via neutrino cooling is not quick enough to allow this material 
to settie quickly so that the shock continues to move outward. 
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Fig. 12. Evolution of tiie central values of the entropy (top) and lepton 
number (bottom) as functions of the central density during the col- 
lapse of Models sl5Nso_ld.b (dash-triple-dotted), sl5Gso_ld.b (dot- 
ted), sl5Gio_ld.b (solid), and sl5Gio_ld.a (dashed). 
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Fig. 13. The top plot shows the shock positions versus time for all pre- 
sented models. Thin lines are for ID models, thick lines for 2D, see 
Fig. 14 for the labelling. The plot below shows the gain radius which 
separates the neutrino cooling from the neutrino heating layer farther 
out. Note that the curves in the lower plot were smoothed over inter- 
vals of 5ms. The evolution of Models sl5Gio_ld.a and sl5Gio_32.a 
was followed to later times. Their complete evolution can be seen in 
Figs. 28, 36, and 42. 



After 6ms, when the shock reaches SOkm, this expansion slows 
down because the mass accretion rate drops and neutrino cool- 
ing has become more efficient, leading to an "equilibrium" be- 
tween mass accretion through the shock and mass settUng onto 
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Fig. 14. The top plot shows the total net heating rate in the gain layer, 
i.e. between the gain radius and the shock. The lower plot shows the 
total cooling rate below the gain radius. Note that the curves were 
smoothed over intervals of 5ms. 
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Fig. 15. The three plots show the neutrinosphere radii of all presented 
models for Ve (top), Ve (middle), and v^_t and v^,^ (bottom), see Fig. 14 
for the labelling. Note that Model sl5Gso_ld.b* is missing in these 
plots (because of data loss). 



the PNS. The shock now slowly expands to 140km at 60ms, 
see Fig. 13. Then, however, the shock slowly retreats again as 
a consequence of the high ram pressure. The neutrino heating 
in a relatively narrow gain layer (see also 5tE^\ in Fig. 14) is 
too weak to support shock expansion to larger radii. At 170ms, 
when the shock has reached a radius of less than SOkm, we see 
a transient reexpansion of the shock, which is stopped again 
after a few km. This feature results from a sudden drop of the 
density in the progenitor (see Fig. 10) and thus of the mass ac- 
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Fig. 16. The three plots show the average energies of the emitted Ve 
(top), Ve (middle), and v^,t or (bottom) as measured by an observer 
at rest at 400km. 
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Fig. 17. Mass accretion rate through the shock for Model sl5Gio_ld.b. 
For the other models its evolution is nearly the same, unless an explo- 
sion occurs. 



cretion rate (Fig. 17) and of the ram pressure, and is connected 
with the "interface" between the silicon shell and the oxygen- 
rich silicon shell (i.e. the Si-SiO interface at IA2Mq). Such a 
phase of transient shock expansion can be seen in all our mod- 
els, with strongly differing consequences (even explosions, see 
Sect. 3.2.2). In Model sl5Gio_ld.b the shock finds a new equi- 
librium radius at only 90km to then follow the subsequent con- 
traction of the nascent neutron star (visible from the decreasing 
neutrinosphere radii, see Fig. 15). 




Fig. 18. Comparison of the neutrino luminosities for all models for an 
observer at rest at 400km. The plot on top shows the total luminosities, 
the three lower ones those for the different v types. The plot at the 
bottom gives the luminosities of v^u, v^, Vt. or v^, which are all assumed 
to be equal. 



The neutrino emission reflects the different phases of the 
evolution: After the collapse phase, during which only Vg are 
emitted in significant numbers, the prompt Vg burst is created 
few ms after bounce (see Fig. 19; note that the signal is delayed 
by approximately 1ms due to the time-of-flight). The burst lasts 
approximately 25ms, has a fwhm of 6ms, and reaches a peak lu- 
minosity of about 350FOE/S (FOE = ten to the fifty-one ergs), 
see Fig. 19. Simultaneously, the newly created PNS starts emit- 
ting Ve and Vx. Their rise in luminosity takes 35ms and 15ms, 
respectively. In the following accretion phase, the emitted neu- 
trinos come from two regions, from a layer below the neutri- 
nosphere "inside" the cooling PNS and to a smaller part from 
the layer of newly accreted material between the PNS surface 
and the gain layer. The luminosity produced by the latter source 
depends sensitively on the mass accretion rate as can be seen 
in the drop of the luminosity when the Si-SiO interface enters 
the cooling region. Interestingly, the mean energies of Ve and Vx 
become very similar after about 100ms post-bounce. This fact 
has been discussed in detail in Keil et al. (2003) and wiU be 
addressed in more detail below (Sect. 3.1.2). 

Looking more closely at our simulation of Model 
sl5Gio_ld.b, we find that its outcome is not very surprising. 
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Fig. 19. The uppermost plot shows the electron neutrino burst, again 
for an observer at rest at 400km. Note the different scale on the ab- 
scissa of this plot. The lower three plots display the neutrino lumi- 
nosities of Ve, Ve, and heavy-lepton neutrinos individually, evaluated at 
their respective neutrinospheres (for a comoving observer). Note that 
the differences compared to the luminosities shown in Fig. 18 come 
mostly from neutrino emission and absorption in the cooling and heat- 
ing layers outside of the neutrinosphere (and not from observer frame 
motion). Also note that Model sl5Gso_ld.b* is missing in the lower 
three plots (due to missing data). 

After shock stagnation and during the whole subsequent evo- 
lution we see that the shock is unable to stop the rapid infall 
of the accreted material: Behind the shock, the matter still has 
negative velocities of several lOOOkm/s, thus falling quickly 
through the narrow gain layer. As an example, at the time of 
maximal shock expansion at 72.5ms, the infall velocity behind 
the shock is still 4 x lQ\m/s (see Fig. 20), and the short dis- 
tance between shock and gain radius of at most 50km corre- 
sponds to an advection timescale through the gain layer of less 
than 15ms. With a moderate neutrino heating rate of 300MeV/s 
per baryon only about 4.5MeV per baryon can be deposited in 
these 15ms, an amount which is clearly insufficient to promote 
shock expansion: At the high entropies (> lOke/by) behind 
the shock nuclei are nearly completely dissociated to nucleons. 
Only a small fraction of a-particles (less than 20%) survive the 
shock passage, and are quickly dissociated during infall. Thus, 
in the shock the infalling matter loses 8-9MeV per baryon due 
to nuclear photo-dissociation. 



3.1 .2. Variations of the input physics 

"Standard" opacities vs. improved opacities We have per- 
formed most of our calculations either using a set of neu- 
trino interactions which we call "standard opacities" (Case so) 
or our fully updated ("state-of-the-art") description of neu- 
trino interactions which we call "improved opacities" (Case 
io)^ The treatment of neutrino interactions in Case io is de- 
scribed in Appendix A. Case so uses the neutrino interactions 
from Bruenn (1985) and Mezzacappa & Bruenn (1993a,b), im- 
plemented as described by Rampp & Janka (2002), but sup- 
plemented by neutrino pair creation and annihilation due to 
nucleon-nucleon bremsstrahlung. In neither of our simulations 
we take into account neutrino pair creation by the plasmon pro- 
cess because this process is negligible for Ve and Ve, and for 

it is either dominated by nucleon-nucleon bremsstrahlung 
(at high densities) or by the flavor coupling vv pair process 
(VxVx <r-¥ VeVe; at low densities). For a recent calculation of 
the plasmon process at SN-relevant conditions see the paper 
by Ratkovic et al. (2003). 

In order to compare the properties of the different inter- 
action rates we have plotted in Figs. 21 and 22 the opacities 
(inverse mean free paths) of all reactions included in Model 
sl5Gio_ld.b for Ve, Ve, and heavy-lepton neutrinos (V;^) for two 
different energies and for the conditions present in the model 
at two different post-bounce times. For completeness, we also 
show the rates for neutrino scattering and absorption on nucle- 
ons as used in Case so, i.e. the rates with the approximations 
of ignoring nucleon recoil and thermal motions, weak mag- 
netism, effective nucleon masses, and nucleon-nucleon corre- 
lations. During the discussion, remember that absorption and 
emission are linked by detailed balance. Phase space blocking 
of neutrinos (and other leptons) in the final channels of the re- 
actions is included in all reactions. 

Charged-current absorption processes with nucleons are the 
dominant mode of destroying Ve and Ve after bounce, when the 
post-shock material is essentially fully dissociated. However, 
we can see that in the PNS core, the improved rates are reduced 
significantly compared to the "standard" rates of Case so, for 
Ve by a factor of up to 30, for Vg even by up to 100. In the 
region around the neutrinosphere, including the gain layer, re- 
coil and weak magnetism reduce the Ve absorption rate by 10- 
20%, while the effects nearly compensate each other in case 
of Ve (compare also Appendix A and Horowitz & Li 1999). 
Absorption on nuclei is dwarfed by nucleon absorption even in 
the pre-shock material. However, our rates for Ve absorption by 
nuclei are still rather approximative (not taking into account re- 
cent advances of electron capture rate calculations by Langanke 
et al. 2001) and Ve capture by nuclei is neglected. Improving 
this treatment is desirable. 

Neutrino-antineutrino pairs can be produced and destroyed 
in different charged- and neutral-current reactions. One of 
these is nucleon-nucleon bremsstrahlung, which contributes 
significantly to the emitted flux of heavy-lepton neutrinos Vx 
(Thompson et al. 2000; Keil et al. 2003; Thompson et al. 2003) 

' The new treatment of electron captures on nuclei during core col- 
lapse with rates from Langanke et al. (2003), however, is not yet in- 
cluded in any of our models presented here. 
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because this process dominates the production of Vx at high 
densities. Interestingly, bremsstrahlung is also the dominant 
absorption reaction for Ve in the PNS core, where the electron 
degeneracy is very high and therefore Ve are orders of magni- 
tudes more abundant than Ve. It is also the dominant produc- 
tion rate since positrons are very rare in the degenerate region, 
leading to a strongly suppressed rate of e^ -h n — > Ve -n p. The 
second pair process is the annihilation of vv pairs to e'^e" pairs. 
The inverse reaction dominates bremsstrahlung in creating Ve 
at low densities and is also a very important process for the 
generation of Vx. This rate, however, is surpassed by the rate 
of vv pair conversion between different flavors, which in most 
regimes is approximately a factor of two more important for 
producing VxVx pairs than e'^e" annihilation. Both rates dom- 
inate the Vx production around the neutrinosphere, thus being 
crucial in the spectrum formation process of these neutrinos 
(cf. Keil et al. 2003; Buras et al. 2003a). For higher energies, 
the leptonic pair production rates become increasingly impor- 
tant relative to bremsstrahlung. 



Finally, scattering processes contribute to the neutrino 
opacity. Scattering of neutrinos on nuclei dominates the scat- 
tering rates above the shock. Scattering on nucleons are the 
dominant rates in the NS core and postshock layer. Dense 
medium effects (see Appendix A), however, reduce the rate 
by a factor of 2-3 at high densities below the neutrinosphere 
and phase space blocking effects lead to a dramatic reduction 
for scattering rates for Ve in case of nucleon recoil being in- 
cluded. Again, weak magnetism and recoil produce a 10% dif- 
ference for Ve-nucleon scattering around the neutrinosphere, 
and a much smaller effect for Ve (see Appendix A). Scattering 
off electrons and positrons and off Ve, Ve are of similar impor- 
tance for Vx (see also Buras et al. 2003a)^. 

^ Note that the strong decrease of the VxV —> v^v rate at large radii 
{i 80km) where the process is essentially unimportant is due to an 
artificial suppression of this rate with a factor of (l + p^^) where 

Pio = p/^lO'" gcm"^^. This procedure was necessary because our 
rate implementation works correctly only when electron neutrinos are 
in chemical equilibrium. 
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Fig. 21. Opacities for the neutrino interactions talcen into account in Model sl5Gio_ld.b at f = 11.6ms post-bounce for two representative 
neutiino energies. Neutrino ("v") reactions with elections ("e", solid), nucleons ("N", dash-dotted), nuclei ("A", dashed), and electron-type 
neutrinos ("Ve", dash-triple-dotted, for only) are shown. For comparison the "standard" opacities for the neutrino-nucleon interactions are 
also displayed (dotted). Thin lines represent scattering processes, thick lines correspond to absorption processes. Note that different from this 
convention nucleon-nucleon bremsstrahlung (important only below the shock) is also represented by thick dashed lines. (See also caption of 
Fig. 22). 



At higher neutrino energies the different opacities show 
similar behaviour (Fig. 21). The suppression of neutrino ab- 
sorption on nucleons at high densities is weaker (because phase 
space blocking is less strong), and bremsstrahlung processes 
are less important compared to other processes. At later times 
the temperature in the postshock layer increases and e" degen- 
eracy is reduced so that the e^e" pair process and ve* scattering 
become more important (compare Fig. 22 with Fig. 21). 



The influence of the neutrino interactions can be 
seen in dynamical simulations with different input physics 
(Figs. 13, 14, 15, 16, 18, 19). In addition to the above presented 
Model sl5Gio_ld.b with our set of "improved" opacities, we 
have performed two dynamical simulations sl5Gso_ld.b and 
sl5Gso_ld.b*, which both were run with the traditional ap- 
proximations for the opacities. In the latter case bremsstrahlung 
was also neglected in order to use the rate input of previous 
simulations by other groups. The influence of nucleon-nucleon 
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Fig. 22. Same as Fig. 21, but for t = 73.4ms after bounce. The opacities include processes with all possible interaction partners: "N" stands 
for nucleons and can represent protons and/or neutrons, "e" can be electrons and/or positrons, "A" represents all nuclei, "v" can be any type of 
neutrino or antineutrino, "Ve" can be electron neutrinos and/or antineutrinos, and "Vx" can be any type of heavy lepton neutrino or antineutrino. 
The opacities include all phase-space blocking factors. The short vertical lines on the top of the plots indicate the radii of the neutrinosphere 
for the given neutrino type and energy (solid) and of the shock (dashed). 



bremsstrahlung was already discussed by Keil et al. (2003) and 

Thompson et al. (2003), the importance of interactions between 
neutrinos of different flavors was pointed out by Buras et al. 
(2003a). 



The collapse phase is hardly influenced by the use of 
our improved description of neutrino-nucleon interactions, be- 
cause free nucleons are not very abundant during core col- 



lapse (Fn + ip - 10"^)^. The improvement of the electron cap- 
ture rate on protons slightly increases the Ve production below 
the trapping density (see Fig. A.4). With the coherent scatter- 
ing rate, which dominates the opacity during infall, not being 
changed, more Ve are created and escape before the matter be- 
comes optically thick. As a consequence, the simulation with 
improved opacities yields slightly lower entropies i\As\ ^ 0.06) 



^ During core collapse until bounce, our simulations were run with- 
out Ve and Vx, both being irrelevant during this phase. 
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Fig. 23. Neutiino flux specti-a for an observer at rest, evaluated at a 
radius of 400km at a time of 30ms after bounce. The diiferent pan- 
els display results for v^, v^, and y>, respectively. The different lines 
correspond to different ID models as labeled in the plot. 



and a slightly decreased lepton fraction (|AFiep| ^ 0.01) in the 
core after trapping compared to the simulation with standard 
opacities (Fig. 12). Therefore the collapse time also decreases 
(|A/coii| - 2.3ms) and the homologous core becomes slightly 
smaller, which reduces the mass enclosed by the shock forma- 
tion radius (lAMsd ^0.02). 

After bounce only the first 50ms are suited for comparing 
the models, in this early phase the structure of the PNS (and the 
above lying layer where the neutrinospheres are located) is still 
very similar in all general relativistic models. 

The lowest panel in Fig. 16 shows that the models cluster 
into two categories concerning the mean Vx energies. The de- 
crease by ~ 2MeV can be explained mainly by the effects of the 
neutrino-nucleon interactions (without bremsstrahlung) when 
energy transfer by nucleon recoil in v-n, p scattering is taken 
into account (see also RafFelt 2001). Bremsstrahlung (com- 
pare models sl5Gso_ld.b* and sl5Gso_ld.b) and the flavor- 
coupling neutrino reactions (compare models'* sl5Gio_ld.a and 
sl5Gio_ld.b) have a small effect on (cy^) (~ O.SMeV). The 
mean electron type neutrino energies are hardly changed by the 
improvements of the neutrino opacities, only the mean energy 
of Ve increases slightly by ~ 0.7MeV. 

Similarly, the Ve spectrum is not altered by the improve- 
ments of the opacities, see Fig. 23. For the high-energy flux 
increases by a factor of three, being consistent with the increase 
of (ey^). The most significant change can again be seen for v^. 
We observe a pinching of the spectrum (higher peak, steeper 
drop at high energies). 



* Note that the omission of jS-dependent terms has no noticeable in- 
fluence on the layers of neutrino emission at early post-bounce times. 



Concerning the luminosities, see Fig. 18, we again find a 
clustering of the models in case of the Vj, luminosity. Here, 
however, the increase by 25-30% comes from the interactions 
between the neutrinos of different flavors, whereas the im- 
proved neutrino-nucleon interactions have a minor (< 10%) in- 
fluence (compare sl5Gso_ld.b* with sl5Gio_ld.a). Again, the 
effects on Vg and Ve are much smaller. For both Ve and Ve the 
luminosities increase by less than 10%. 

As a direct consequence of the increased luminosities, the 
cooling rate of the PNS, dtE^ooX in Fig. 14, is 25% higher in 
case the interactions between neutrinos of different flavors are 
turned on. The effect on the evolution of the shock, however, 
is very weak and becomes visible only later than 100ms after 
bounce when the enhanced contraction of the PNS in models 
including the flavor-coupUng neutrino interactions has altered 
the structure of the star between neutrinosphere and shock. 

Newtonian vs. general relativistic simulations Finally we 
discuss the differences arising between simulations with 
our approximative treatment of general relativity (GR) and 
Newtonian calculations by comparing Models sl5Gso_ld.b 
with sl5Nso_ld.b. In the Newtonian simulation, the gravita- 
tional potential is less deep and the onset of core collapse is 
not supported by relativistic destabilization. Hence the collapse 
time (until shock formation) of 196.1ms is longer compared to 
177.7ms for the relativistic potential. Nevertheless the central 
lepton fraction after trapping is nearly the same (it is smaller 
by only AFiep - 0.002 in the Newtonian case, see Fig. 12), 
suggesting a very similar enclosed mass at shock formation. 
However, due to the lower gravitational forces and infall veloc- 
ities the shock forms at an enclosed mass of Mgc = O.62M0 
and a radius of = 12.5km in the Newtonian case, whereas 
it is launched at Mjc = 0.5 IM© and r^^ = 10.8km in the GR 
simulation. 

The less compact structure of the PNS in the run with the 
Newtonian gravitational potential strongly influences the post- 
bounce evolution: The neutrinospheres and the accretion shock 
stand at larger radii, the increase of temperature and entropy in 
the gain layer is slower, and the neutrino luminosities are lower, 
see Figs. 13, 15, 18, and 24. These results agree with the more 
elaborate discussion of differences between Newtonian and GR 
simulations by Liebendorfer et al. (2001). 

3.1 .3. Velocity-dependent terms in the first order 
moments equation 

Terms which contain /? in the first order moments equation of 
the transport equation (except for the ^d/dr-terms) are usu- 
ally considered to be small and therefore negligible for many 
problems when v/c «: 1 (Mihalas & Mihalas 1984; see also 
the discussion in Rampp & Janka 2002). However, in the pre- 
sented supernova simulations, we found that the velocities are 
sufliciently large in the neutrino decoupling and heating lay- 
ers, which are crucial for the explosion, that these "j8-terms" 
can make a big effect. With "y6-terms" we mean here all terms 
depending on /3r, except the advection terms, which are pro- 
portional to fird/dr and which are included in all calculations. 
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Fig. 24. Comparison of radial profiles for different variables at post-bounce times of 12.3ms (solid) and 177.7ms (dashed) for Models 
sl5Nso_ld.b (thin) and sl5Gso_ld.b (thick). L,ot is the total neutrino luminosity, i.e. the sum of contributions from neutrinos and antineu- 
trinos of all flavors. 



The omission of these yS-terms in the neutrino momentum equa- 
tion turned out to change the dynamical evolution of the su- 
pernova and to lead to an artificial explosion in ID and 2D 
(Sects. 3.1.4 and 3.2.2). Having in mind that in the Oiyjc) 
approach all velocity-dependent first-order terms in the radi- 
ation energy equation are known to be relevant (see Mihalas 
& Mihalas 1984) and the importance of the j6^-terms was rec- 
ognized in Sect. 2.3.2, we therefore conclude that all terms 
of first order in the fluid velocity must be taken into account 
in supernova simulations^. Doing so, we found that our trans- 
port code produces results in very good agreement with spheri- 
cally symmetric simulations with the Oak Ridge/Basel AGILE- 
BOLTZTRAN code for Newtonian as well as relativistic grav- 
ity (Liebendorfer et al. 2005; Marek et al. 2005). 

We want to investigate the role of the /?-terms in the neu- 
trino momentum equation, Eq. (B.14), more closely with a 
snapshot of Model sl5Gio_ld.b at t^\, = 114ms, for which we 
have evaluated these terms at energies near the spectral maxima 
of neutrinos and antineutrinos, see Fig. 25. For Ve and Ve the yS- 
terms account for up to 20-30% of the LHS of Eq. (B. 14) in the 
neutrino-heated region between 70 and 1 10km; for muon neu- 
trinos, which do not contribute to the heating at a significant 
level, the jS-terms can even dominate the LHS of Eq. (B.14). 
The individual size of these terms at 1 14 ms after core bounce 
is visible in Fig. 25d. One can see that the different terms have 

' We point out that some authors think that going even beyond the 
Oiv/c) approach of our work is desirable (cf. Cardall et al. 2005). 



different relative importance in different regions of the star. We 
emphasize that it can therefore not be concluded that any of 
these terms can be safely ignored for all situations. 

The role of the j8-terms can be made more transparent if we 
write them on the RHS and compare them with the source term 
for momentum exchange. In regions with negative velocity and 
negative velocity gradient, the yS-terms have the opposite sign of 
the source terms in Eq. (B.14) and therefore effectively reduce 
the coupling between matter and neutrino flux. Thus, neglect- 
ing the j8-terms has the opposite effect, leading to a less rapid 
streaming of neutrinos and an increased energy density of neu- 
trinos above and below the gain radius. Note that for higher 
energies, the velocity terms become less important. 

These statements are confirmed by Fig. 26, where we show 
stationary transport solutions calculated for the considered stel- 
lar profile both with and without the j0- terms. The neutrino flux 
H and the flux factor HjJ are larger in the solution with the 
jS-terms, and the neutrino density is smaller, which directly de- 
creases the neutrino heating (Fig. 26b). From the plots, it can 
also be seen that the solution with the yS-terms is consistent with 
physical requirements: Since the neutrino quantities are given 
in a comoving frame, both luminosity and energy density re- 
veal a larger blue- shifting ahead of the shock where the infall 
velocities are higher. The corresponding "step" at the shock, 
however, is absent in case of J in the solution without y6-terms 
(Fig. 26a, left panel). The flux factor is affected by these differ- 
ences. Without yS-terms it clearly disagrees with results of trans- 
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Fig. 25. Comparison of the different terms in the neutrino momentum equation (Eq. B.14, multiplied by mty/p) for Model sl5Gio_ld.b at 
= 114ms for Ve (top left), Vc (top right), and heavy-lepton neutrinos, (bottom left), at the indicated energies. The different lines represent 
the j0-terms (thick solid), the remaining terms on the LHS of Eq. (B.14) (dotted), and the source term on the RHS (dashed). We have also 
plotted the sum of all three (dash-dotted), which reveals minor deviations from zero as a consequence of the post-processing evaluation of 
the discretized gradients. Finally, the thin solid line in the upper left plots shows the velocity profile. Panel d displays for Vg the different j8r- 
terms in the neutrino momentum equation (Eq. B.14, multiplied by wiby/p) individually. These correspond to the physical effects of radiation 
compression, Hr^^drir^jir) (thick dashed) and the term HdrPr (thin dashed), Doppler effects (thick dash-dotted), and fluid acceleration (thin 
dotted). 



port calculations with an >S;v solver for the Boltzmann equa- 
tion (see Liebendorfer et al. 2005). In contrast, the variable 
Eddington factor /ed = KIJ turns out not to reveal the same 
degree of sensitivity to the yS-terms. 

The yS-terms were found to have a dramatic dynamical ef- 
fect. In Fig. 13 we see the evolution of Models sl5Gio_ld.a 
and sl5Gio_ld.b, which mostly differ in the fact that only 



Model "b" includes the yS-terms*'. We see that the shock po- 
sitions evolve increasingly differently. Initially, in Model "a" 
the post-shock heating by neutrinos is higher, see Fig. 14. At 
later times, the increased heating has driven the shock farther 
out, thus increasing the size and mass of the gain layer. The 



* The disregard of flavor-coupling neutrino processes in Model 
sl5Gio_ld.a is not very relevant in the present context. 




Fig. 26. a Comparison of steady-state transport solutions without (left) and with the j8-terms in the neutrino momentum equation, computed 
for a background profile from Model sl5Gio_ld.b at tpb = 114ms. The neutrino energy density J (times r^), the flux factor = H/J (both 
with thick lines), the neutrino energy flux H (times r^), and the variable Eddington factor /ej = K/J (both thin) are given for electron neutrinos 
(solid), electron antineutrinos (dashed) and muon neutrinos (dash-dotted) for an observer comoving with the stellar fluid. The kink of fu within 
the shock is a numerical artifact produced by the interpolation of J and H between the grid zones of a staggered mesh, b Corresponding 
velocity profile (thick), total neutrino luminosities (as the sum of the luminosities of neutrinos of all flavors; dashed) and the energy source 
terms (neutrino heating minus cooling; thin solid Une). 



resulting positive feedback between increasing gain layer and 
thus increasing heating on the one hand and expanding shock 
radius on the other leads to a much larger shock radius than in 
Model "b", where the contracting shock causes the gain layer 
to become more and more narrow. During the phase of rapid 
accretion until about 180ms post-bounce, the total heating rate 



in the gain layer is larger in Model sl5Gio_ld.b by a factor of 
about two, and even larger at later times. 

This analysis is confirmed when Models sl5Gio_ld.b (with 
y6-terms) and sl5Gio_ld.a are compared at f = 73ms post- 
bounce (Fig. 27): The flux factor = H/J in the former model 
increases much faster in the regions of effective neutrino cool- 
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Fig. 27. a Transport results at tpb = 73ms for Model sl5Gio_ld.a where the j8-terms were dropped, and for Model sl5Gio_ld.b, which was 
computed with all terms in the neutrino momentum equation. The neutrino energy density J (times r^) and the flux factor fu = H/J are plotted 
with thick lines, the neutrino energy flux H (times r^) and the variable Eddington factor = K/J with thin lines for electron neutrinos (solid), 
electron antineutrinos (dashed) and muon neutrinos (dash-dotted), b The corresponding profiles of velocity (thick), total luminosity of neutrinos 
and antineutrinos of all flavors (dashed) and the source term for energy exchange between neutrinos and stellar medium (thin solid). 



ing and heating than in Model sl5Gio_ld.a. The lower neutrino 

energy density corresponds to a significantly reduced integral 
rate of neutrino heating and a much reduced rate of cooling 
below the gain radius. 



3.1 .4. A model with marginal behaviour 



Although Model sl5Gio_ld.a was computed with reduced ac- 
curacy in the neutrino transport (the yS-terms in the first order 
moments equation were neglected), its evolution deserves a dis- 
cussion. It reveals an oscillating shock behaviour with growing 
amplitude (similar to what was previously seen in models of 
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Fig. 28.Post-bounce evolution of Model 
sl5Gio_ld.a. The upper panel shows the 
summed comoving-framc luminosity of 
electron neutrinos and antineutrinos at the 
neutrinosphere of Ve (soUd) and at the 
gain radius (dotted), and this luminosity at 
400km for an observer at rest (dashed). The 
lower dashed line corresponds to the (indi- 
vidual) luminosity for neutrinos of p or t 
type, also for an observer at rest at 400km. 
The middle panel displays the evolution of 
the shock (upper thick line), gain radius, 
and Ve-sphere (lower thick line). The thin 
Unes represent mass shells, starting at M = 
1.41Mo and increasing in mass in steps of 
0.002Mo. The lines were truncated at the 
Ve-sphere. The lower panel depicts the to- 
tal energy loss rate in neutrinos from the 
cooling layer (solid, divided by 10) and the 
net energy deposition rate by neutrinos in 
the gain layer (dashed). Note that here, in 
contrast to Fig. 14, the curves were not 
smoothed over time intervals. 



the Livermore group, see e.g. Wilson et al. 1986), which finally 
leads to an explosion. 

As discussed in Sect. 3.1.3, the omission of the yS-terms 
causes increased neutrino heating behind the stalled SN shock. 
The shock is therefore driven farther out before it stagnates 
at about 180km. When the interface between silicon layer and 
oxygen-rich silicon shell arrives at the shock, a sudden decrease 
of the ram pressure occurs and instigates a sequence of peri- 
odic phases of shock expansion and contraction (Figs. 28, 29). 
Whenever the shock expands and the mass infall is decelerated, 
the neutrino luminosity drops with a slight delay as a conse- 
quence of the decreased mass accretion rate towards the cool- 
ing region. In response to the decreased luminosity the heating 
almost simultaneously drops, partly also because the accreted 
material expands behind the shock and stays away from the 
gain radius where heating is most efficient. With the energy in- 
put being quenched the matter behind the shock quickly falls 
inward. Being drained of its pressure support the shock expan- 



sion comes to a halt again and the shock retreats. The sudden 
increase of the mass flow into the coohng region raises the ac- 
cretion luminosity extremely, and the heating rate reaches up to 
20% of the cooling rate, also enhanced by more mass moving 
into the region just above the gain radius. The shock bounces 
and is driven out in a new phase of expansion. This cycle is 
repeated several times. Obviously, the feedback of the cycle is 
positive, since the ampUtude of the shock expansion increases 
with time and the cycle finally leads to an explosion. 

It is beyond the scope of this work to determine exactly 
the criterion when such a feedback cycle is obtained and why 
Model sl5Gio_ld.a but no other of our models shows such a 
vertiginous fate. The fact that the yS-terms are neglected cer- 
tainly increases the heating, and a sufficiently large driving 
force seems crucial for entering the oscillatory mode. 

We see two possibilities for explaining the phenomenon, 
however make no effort to perform a detailed analysis of any of 
these. On the one hand the physical conditions which charac- 
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Fig. 29. Radial structure of Model sl5Gio_ld.a at post-bounce times of 228.6ms (solid), 236.4ms (thin dashed), 254.4ms (thick solid), and 
262.6ms (thick dashed). At the first and last time, the post-shock velocities reach their maximum positive values, at the second time the shock 
reaches its local maximum, and at the third time the post-shock velocities reach their minimum negative values. At 254ms the maximum cooUng 
rate reaches Qe = -30.1GeV/s per baryon. 



terize the expansion and contraction phases resemble the con- 
ditions that allow for non-adiabatic vibrational instability (k- 
mechanism) in stellar atmospheres where the outward-going 
heat flow is modulated by the rhythm of the pulsation and in- 
stability occurs if during adiabatic compression the absorption 
coefficient increases (see, e.g. Kippenhahn & Weigert 1990, 
Chapter 39.1). 

On the other hand, the model might reveal the action of 
the so-called acoustic-advective cycle proposed by Foglizzo & 
Tagger (2000) for adiabatic accretion flows to black holes. In 
this scenario, acoustic waves created at the PNS surface prop- 
agate to the shock and cause entropy fluctuations there. When 
these fluctuations are advected to the PNS surface by the ac- 
cretion flow, they create new acoustic waves and so on. If the 
feedback is positive, an / = mode can be built up, as seen in 
our model. A preliminary analysis by T. Foglizzo (2003, per- 
sonal communication) revealed that our oscillating model has 
favorable conditions for developing an / = instabiUty by the 



advective acoustic cycle. However, in our model neutrino heat- 
ing and cooling dominate and the accretion flow is not at all 
adiabatic. A more reliable analysis has to account for this fact 
and remains to be done. 

3.2. Two-dimensional models 

An inspection of the stability criterion (see Sect. 2.3.2 for our 
specific definition of the stabiUty criterion) for the ID mod- 
els (Fig. 30) tells us that the gain layer as well as an ex- 
tended region in the PNS are convectively unstable, so that 
multi-dimensional simulations become mandatory. Here, we 
present the first two-dimensional simulations of convection 
in supernova cores using a multi-frequency treatment of neu- 
trino transport. The two Models sl5Gio_32.a and sl5Gio_32.b 
correspond to the one-dimensional Models sl5Gio_ld.a and 
sl5Gio_ld.b, respectively, presented in Sects. 3.1.1 and 3.1.4. 
Both 2D simulations include our most elaborate descrip- 
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tion of neutrino-matter interactions (Case "io") as described 
in Appendix A of this paper and an approximative, spher- 
ically symmetric treatment of general relativity (Sect. 2.1). 
Remember that Model "b" is the one with the full treatment 
of neutrino transport, while Model "a" incorrectly neglects 
the jS-terms in radial direction (see Sect. 3.1.3). In the angu- 
lar direction both models employ 32 zones with a resolution 
of 2.7°, covering a wedge around the equator from -43.2° to 
-1-43.2°. Both models were followed through core collapse in 
one dimension and mapped to the 2D grid around 7ms after 
the bounce. Simultaneously, perturbations were seeded by ran- 
domly changing the radial velocity by up to +1% in the entire 
star. The treatment of the lateral neutrino transport was carried 
out as described in Sect. 2.2. Note that in Model sl5Gio_32.b 
the inner 2km were treated spherically symmetrically to main- 
tain a reasonably large CFL time step. Model sl5Gio_32.a did 
not include the lateral momentum transfer to the fluid, Eq. (18), 
for which reason the core inside a radius of 25km had to be cal- 
culated in spherical symmetry. Thus artificial hydrodynamic 
instability could be avoided, but also the physical PNS con- 
vection would not be followed in a layer below the neutri- 
nospheres. Note that the entropy wiggles in the postshock layer 
discussed in detail in Sect. 2.5 were present in both 2D simu- 
lations presented here and introduce artificial perturbations of 
the stellar models in radial directions. More recent calculations 
with our improved scheme of treating rest-mass effects in the 
EoS (Sect. 2.5 and Appendix C), however, showed that these 
numerical perturbations on the entropy profile did not have an 
appreciable influence on the development and growth rates of 
convection in the hot bubble region. 

3.2.1. s15Gio_32.b: A model with full transport 
treatment 

Model sl5Gio_32.b represents a 2D supernova computation 
with the presently most complete and detailed implementation 
of spectral neutrino transport in such models. It fails to explode. 
In Fig. 31 we show the evolution of shock trajectory and "mass 
shells" of this model. Note that a "mass shell" trajectory of a 
2D simulation does not represent the evolution of Lagrangian 
mass elements as in spherical symmetry but describes the time 
evolution of the radius of the spherical shell enclosing a certain 
value of the mass. After the prompt expansion, the shock is 
pushed slowly outward by accumulating mass until it reaches a 
maximum radius of 150km around 60-70ms, and then starts 
to slowly retreat again. This is quite similar to the ID case 
(c.f. Fig. 13). Convection in the neutrino-heated layer and in 
the PNS shifts the shock position farther out by just several 
10km. 

The initial hope in performing 2D calculations was that 
the "hot bubble" (HB) convection in the gain layer would 
strengthen the shock sufliciently to lead to an explosion, sim- 
ilar to what was found previously by Herant et al. (1994) 
and Burrows et al. (1995) with neutrino transport being de- 
scribed by grey, flux-limited difiiision. Hot bubble convection 
was found to be helpful for shock revival by means of trans- 
porting energy to the shock, thus enhancing the postshock pres- 




Fig.30. a Brunt- Vaisala frequency as defined in Sect. 2.3.2 at the 
post-bounce times of 72.3ms (dashed) and 181.4ms (soUd) for Model 
slSGio.ld.b (thin) and sl5Gio.32.b (thick). The lines are plotted only 
downstream of the shock. Positive values indicate instability. For the 
2D model the evaluation was performed with laterally averaged quan- 
tities, b Same as a, but for models sl5Gio_ld.a (thin) and sl5Gio_32.a 
(thick) at post-bounce times of 72.2ms (dashed) and 150.1ms (soUd). 

sure and increasing the efficiency of neutrino energy deposi- 
tion: Strongly heated material close to the neutrinosphere is 
carried to the shock in bubbles rising due to buoyancy forces. 
At the same time narrow downflows passing besides the bub- 
bles feed the region of strong neutrino heating with cool mate- 
rial. 2D calculations also get rid of another big disadvantage of 
ID calculations: When the shock starts expanding in ID calcu- 
lations, the accretion of fresh material into the cooling region 
is quenched, the sudden drop in luminosity and thus heating re- 
duces further support for the expanding shock and thus can be 
fatal for an explosion. 

In the current model, although the whole gain layer is con- 
vectively unstable from fpb = 40ms on (see Fig. 32a), the con- 
vection has difficulties developing: Becoming distinct at 60ms 
after bounce, the HB convection (Fig. 33) has not yet grown 
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Fig.31."Mass shell" trajectories for Model 
sl5Gio.32.b. The shells plotted with bold solid 
lines are labelled by the corresponding enclosed 
masses. The plot also shows the transport neutri- 
nospheres for Ve (thick solid), Ve (thick dashed), 
and Vx (thick dash-dotted), the mass shell at which 
the silicon shell becomes oxygen-rich (knotted 
solid line, at 1.42Mq), and the angle-averaged 
shock radius (solid line with dashes). Further, we 
have marked the regions which contain a mass 
fraction of more than 60% in Fe-group elements 
(dark shaded) and of more than 30% in ''He (light 
shaded, visible in the vicinity of the shock). Finally, 
the inner thin dashed line marks the gain radius, 
while the outer one marks the interface between the 
regimes of the low-density and high-density EoSs 
(p = 6 X 10^ gem"'). 



strong and formed large-scale structures when the shock be- 
gins retreating again. Therefore HB convection is unable to 
push the shock significantly farther out relative to the shock 
in the ID simulation. Also during the later evolution the HB 
convection cannot develop strength because within only a few 
milliseconds the matter accreted by the shock is advected into 
the cooling layer and buoyancy forces are not strong enough 
to enforce significant rise against the rapid infall of the post- 
shock flow. The shock itself shows almost no response to the 
convection and remains spherically symmetric throughout the 
run (Fig. 33). The HB convection in this model leads to only 
small quantitative difl'erences compared to the one-dimensional 
results. 

Eventually, the transient shock expansion correlated with 
the sudden decrease of the ram pressure when the oxygen-rich 
silicon shell meets the shock allows the hot bubble convection 
to strengthen because the accretion velocities behind the shock 
are reduced. At this time (180ms), however, the shock has al- 
ready retreated to below 100km, and the infall velocity is too 
large to be inverted by the outward acceleration of buoyant, 
neutrino-heated matter. 

Fig. 32b is suitable for showing the effects arising from 
convection. While entropy fluctuations develop due to the con- 
vective activity after convection has become distinct about 
60ms after bounce, convection is too weak to noticeably per- 
turb the entropy profile about 80ms later. The perturbations 
grow again when convection is transiently revived after the Si- 
SiO interface has passed the shock. Interestingly, the entropy 
perturbations persist for some time when they are advected 
downstream into the convectively stable cooling region (com- 
pare with Fig. 32a). 

A feature visible in Fig. 32 betwen 30 and 50km needs 
to be explained: Large lateral velocities associated with mat- 
ter advected downward from the HB to a convective zone in 
the PNS show up between 100 and 150ms post-bounce. This 
phenomenon is caused by the periodic lateral boundary condi- 
tions applied in our simulation, which allow rings of uniform, 



lateral velocity to develop and to be stable. In case of reflect- 
ing boundary conditions such an eff'ect would be suppressed. 
This reminds us that in two-dimensional simulations the cho- 
sen lateral boundary conditions determine the results to some 
extent. In the present case the kinetic energy stored in the rings 
of lateral motion is negligible compared to the energy needed 
to trigger an explosion and we are therefore confident that this 
numerical phenomenon has no macroscopic dynamical conse- 
quences. 

Fig. 32 shows that a second region of convective activity 
is present inside the PNS and persistently encompasses a layer 
between 10-15km and about 30km. The layer of large lateral 
velocities is somewhat more extended in the radial direction 
than the layer of convective instability due to over- and un- 
dershooting of buoyant matter. The absence of visible entropy 
fluctuations in Fig. 32b is explained by the very high efficiency 
with which the PNS convection transports and redistributes en- 
tropy. Therefore the entropy gradient becomes extremely flat in 
this region. Despite of leading to changes of the PNS structure 
and affecting the neutrino emission, this convection inside the 
PNS has no major consequences for the NS evolution during 
the simulated periods of time. We will discuss these effects in 
greater detail in our forthcoming Paper II. 

3.2.2. s15Gio_32.a: A model with artificial explosion 

A comparison of the simulations performed for this work 
shows that none of our models with the full implementation of 
relevant physics (state-of-the-art description of neutrino-matter 
interactions, general relativistic gravity, variable Eddington 
factor closure of coupled set of neutrino moments equations 
and model Boltzmann equation) yields explosions, neither in 
spherical symmetry nor in 2D with convective processes taken 
into account. But as we already have seen in connection with 
Model sl5Gio_ld.a, the omission of the /?,.-dependent terms 
(but not the advection term jSrdjdr, which is always included) 
in the neutrino momentum equation (Eq. B. 14) makes a big dif- 
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Fig. 32. a Convective regions for Model sl5Gio_32.b. Dark shad- 
ing indicates regions where the lateral velocities reach more than 
700km/s, light shading regions where the stability criterion described 
in Sect. 2.3.2 predicts instability. The contour lines correspond to the 
boundaries of the dark-shaded regions. The evaluation was performed 
only downstream of the shock (upper thick solid line) with polar aver- 
ages of the variables. The lower thick solid line depicts the gain radius. 
The thin lines depict the neutrinospheres for (solid), Ve (dashed), 
and Vv (dash-dotted). Finally, the dotted lines show where the densities 
are 10'^ gcm"^ (upper line) and 10'^ gcm"^ (lower line), b Same as a 
but with the light shading marking regions where the standard devia- 
tion of entropy variations (defined in analogy to Eq. 24) is cTj > 0.02. 



ference and favors explosions. Model sl5Gio_32.a is the coun- 
terpart of this ID model in two spatial dimensions and accounts 
for the fact that the neutrino-heated "hot bubble" region behind 
the shock is convectively unstable according to the Ledoux cri- 
terion (Fig. 30b). However the reader should keep in mind that 
the omission of the yS-terms in Eq. (B.14) is an approxima- 
tion which alters the neutrino transport in an unacceptable way. 
Thus, the explosion obtained in Model sl5Gio_32.a has to be 
considered as artificial. 

Note that Model sl5Gio_32.a was calculated with a 2D po- 
lar grid for r > 25km, but spherical symmetry was assumed for 



10 
8 
6 
4 
2 



70ms ^^^^150ms 



14 




12 




10 








8 


\ 


6 




4 


in 


2 










140 



100 60 
r [km] 



20 



40 80 
r [km] 




Fig. 33. Snapshots of Model sl5Gio_32.b with convection in the "hot 
bubble" region. The white lines mark the shock position (thick), the 
gain radius (outer thin), and the Vj-sphere (inner thin). The solid black 
line denotes the equatorial plane of the polar grid used in the simu- 
lation, a These panels show colour coded the entropy s. b All panels 
show colour coded the radial velocity v, . 



radii smaller than 25km, thus ignoring the existence of a con- 
vectively unstable layer below the neutrinosphere. However, 
convection inside the PNS does not have a significant influ- 
ence on the explosion which develops in the discussed model; 
it is important for the cooling of the nascent PNS only on a 
timescale much longer than 100ms. 

For continuing the simulation to late stages with accept- 
able requirements of CPU time we have mapped the 2D cal- 
culation on a ID grid after 468ms of post-bounce evolution. 
At this time, multi-dimensional processes (transient downflows 
of low-entropy matter) affect the onset of the neutrino-driven 
wind below 200km, while the immediate postshock layer is 
convectively only weakly unstable (Fig. 34a) and convective 
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Fig. 34. a Convective regions for Model sl5Gio_32.a. Dark shading 
indicates where the mass inflow rate (evaluated for zones where v,- < 
0) exceeds 10""Mo/s. Light shading denotes regions where aisy > 
30/s (see Sect. 2.3.2 for the definition of the stability criterion). The 
contour lines mark the boundaries of the dark-shaded regions. The 
evaluation was performed only downstream of the shock (upper thick 
solid line) with polar averages of the variables. The lower thick solid 
line corresponds to the gain radius. Note that the core of this model 
inside of 25km was computed in spherical symmetry, b The solid line 
represents the mass outflow rate (evaluated for zones where > 0) 
at r = 100km versus time. The dashed line shows the corresponding 
mass inflow rate. The lines were smoothed over time intervals of 5ms. 



structures in the bulk of expanding matter between this layer 
and the developing wind are going to freeze out and expand es- 
sentially self-similarly. The mapping required to set the lateral 
velocity components as well as the corresponding momentum 
and kinetic energy terms to zero. Along with this change we 
turned on the y6-terms in the neutrino momentum equation to 
follow the evolution of the neutrino-driven wind from the proto 
neutron star more realistically, and switched to the more so- 
phisticated treatment of nuclear statistical equilibrium at low 
densities by a 4-species Saha solver as described in Sect. 2.4. 
The latter replacement was necessary because the previous 




Fig. 35. 2D snapshots of Model sl5Gio_32.a for four different post- 
bounce times, a shows the entropy s, b shows the electron fraction 
Y^. The diagonal black lines mark the equator of the polar grid, the 
horizontal axis corresponds to cos(i?) = -0.67, the perpendicular axis 
to cos()?) = +0.67. The shock radius is indicated by a thick black or 
white line, and the gain radius and Ve-sphere by thin white lines. 



approximative description of NSE in the low-density EoS of 
Rampp & Janka (2002, Appendix B) did not solve the Saha 
equations but described nuclear dissociation and recombina- 
tion by temperature-dependent composition changes which do 
not guarantee entropy conservation and thus cause numerical 
artifacts during the later phases of wind evolution. 

Note that the steepening density gradient at the PNS sur- 
face with ongoing time made it necessary to refine our Eulerian 
grid in this region. Rezoning was done at 206ms, 468ms, and 
759ms. We should also mention that Model sl5Gio_32.a did 
only take into account Si burning, but not C, O, Mg, and 
Ne burning as described in Rampp & Janka (2002). A post- 
processing estimate, however, showed that oxygen burning 
would have contributed less than 0.011 X lO^^erg to the explo- 
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Fig. 36. "Mass shell" trajectories for the explod- 
ing 2D Model sl5Gio_32.a. "Shells" plotted by 
bold lines are labelled with their corresponding en- 
closed masses. Also shown are the transport neutri- 
nospheres for (thick solid), (thick dashed), and 
Vx (thick dash-dotted), the interface between the sil- 
icon shell and the oxygen-rich silicon shell (knotted 
solid line, at 1.42Mq), and the shock and reverse 
shock or wind termination shock (thick solid lines 
with dashes). Moreover, the regions are marked 
which contain a mass fraction of more than 60% of 
Fe-group elements (dark shaded) or more than 60% 
of '*He (shaded). Light shading indicates other re- 
gions which have a mass fraction of more than 30% 
of '^He. The lower dashed line denotes the gain ra- 
dius and the upper one marks the interface between 
the low-density and high-density EoSs (i.e. this Une 
corresponds to p = 6 x 10' g cm"^). 




Fig. 37. Time evolution of the total mass in the gain layer and the 
average local specific binding energy in the gain layer for four models, 
including the exploding Model sl5Gio_32.a. 
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Fig. 38. Explosion energy for Model sl5Gio_32.a as defined in 
Eq. (29). The small jump at fpi, = 468 ms is a consequence of the map- 
ping of the 2D model to a ID grid, and the slow increase afterwards is 
associated with the power of the neutrino-driven baryonic wind. 



sion energy, while the ignition conditions for the other burning 
processes were not reached in our model. 

Evolution Model sl5Gio_32.a develops powerful convection. 
Like in Model sl5Gio_32.b, convective activity in the gain 
layer starts at about 60ms after bounce (see Fig. 35). However, 
in this model neutrino heating is stronger (the heating integral 
is higher roughly by a factor of 2, cf. dfE^i in Fig. 14) due to 
the oinission of the yS-terms, leading to a steeper entropy gradi- 
ent, which in turn means a stronger destabilization of the gain 
layer. Furthermore, the stronger heating is sufficient to prevent 
the shock contraction seen in Model sl5Gio_32.b after 70ms 
post-bounce (see Fig. 14). As a consequence the convective re- 
gion between gain radius and shock remains extended and the 
convective cells can grow and merge to larger structures, lead- 



ing to the largest scale possible on the grid for the boundary 
conditions imposed in our simulation (i.e. a wedge of +43.2° 
around the equator with periodic boundary conditions): After 
about 150ms postbounce one big convective cell with a rising 
bubble and one downflow has formed and transports energy 
efficiently to the shock (Fig. 35). Convection in this model be- 
comes strong enough to push the shock farther out; the gain 
layer, and therefore the convective region, continuously expand 
after convection has set in (Figs. 13, 36). Finally, convection 
supports an explosion, although the corresponding ID model, 
sl5Gio_ld.a, did not explode at such an early time. 

The onset of the explosion occurs when the infalUng 
oxygen-rich silicon shell arrives at the shock at 260km and 
150ms after bounce. The interface between the silicon layer 
and the O-rich sihcon shell is characterized by a steep entropy 
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jump and the preshock density drops to a much lower value. 
This leads to a sudden decrease of the mass accretion rate and 
thus of the ram pressure. The shock reacts to that by expansion 
and rapid acceleration. Although the explosion starts at this 
moment, the existence of the composition interface may not 
be crucial for the explosion of Model sl5Gio_32.a. Favorable 
conditions for an explosion can be seen even before the com- 
position interface has reached the shock: For one the shock is 
already expanding (Fig. 36). Second, the mass contained in 
the gain layer remains constant between 70 and 170ms after 
bounce in contrast to all other simulations, including Model 
sl5Gio_ld.a, where Mgi is always decreasing, see Fig. 37. 
(Note that the difference in Mgi between Models sl5Gio_32.a 
and sl5Gio_ld.a during the first 70ms after bounce arises from 
the omission of the w <-> vv rates in the latter model). Third, 
the matter in the gain layer is accumulating energy: The whole 
gain layer is continuously getting closer to become gravitation- 
ally unbound (Fig. 37) and the energy contained in grid zones 
with positive local specific binding energy ("explosion en- 
ergy"; Eq. 29) already starts growing after 120ms post-bounce 
(Fig. 38). Also, looking at the profiles at /pb - 150ms in Fig. 39, 
one can see that the jump in entropy and density has not yet ar- 
rived at the shock, but that the average radial velocities behind 
the shock have almost become positive. All these facts suggest 
that the convective overturn in the gain layer has been success- 
ful in transporting energy from the gain radius to the shock and 
thus has increased the heating efficiency. 

Shortly after the launch of the outgoing SN shock the den- 
sity in the region between dense post-shock shell and PNS 
drops (Fig. 40) and the gain radius retreats to a position just 
outside the neutrinosphere, which it reaches about 300ms after 
bounce (Fig. 36). The onset of the neutrino-driven wind after 
the start of the explosion is visible in Fig. 40 as a rapid (in 
our simulation spherically symmetric), baryonic outflow with 
velocities up to ~ 5 x lO^cm/s. It develops after ~ 500ms 
post-bounce (see also Fig. 36) and catches up with the slower 
and denser shell of early SN ejecta. In the deceleration region 
a reverse shock (wind termination shock) develops (Janka & 
Miiller 1995; Tomas et al. 2004) which is a very prominent fea- 
ture in the density, velocity, temperature, and entropy profiles 
of Fig. 40 after ~ 600ms post-bounce (see also Fig. 36). 

Figure 34a shows the convective activity during the 2D 
phase of the simulation. While convection in the PNS is in- 
hibited to develop despite of convective instability because 
we calculate with a spherically symmetric grid below 25km, 
large parts of the gain layer are unstable and show convec- 
tive motions (Fig. 35). On the other hand we see that in the 
matter swept up by the expanding shock convective instability 
is so weak (the condition wbv < 30/s appUes) that convec- 
tive activity between about 1000 and 3000km does not set in. 
Around 320ms after bounce the downflows of mass to the PNS 
cease (see also Fig. 34b) and the neutrino-driven baryonic out- 
flow begins to develop (Fig. 36). However, a massive fallback 
in the form of several downflows occurs between 390ms and 
440ms after bounce and destroys this neutrino-driven "wind". 
The wind begins to recover again afterwards (Fig. 36). In 
Fig. 34b we see that the mass inflow (at 100km) becomes al- 
most equal to the mass outflow after the onset of explosion. 



meaning that the net mass flow is small compared to the over- 
turn. Consistent with Fig. 34a, the outflow dominates between 
300ms and 390ms post-bounce and after 440ms post-bounce, 
while the transient downflows between 390ms and 440ms lead 
to strong overturn activity also in the vicinity of the nascent 
NS. 

Explosion energy To calculate the explosion energy at a given 

time, we add up the local specific binding (LSB) energy ej*^^', 
defined in Eq. (27), for all grid zones between neutrinosphere 
and shock where this energy and the radial velocity are both 
positive: 

^expi= Yj &';i»)p{r,mVir,&). (29) 

r,tf,cond 

Note that the gravitational potential cD^™'. enclosed ^^^^ 
calculate the LSB energy was calculated assuming a spheri- 
cally symmetric mass distribution (2D efl'ects are only caused 
by convective mass motions and the associated density fluctu- 
ations are modest and therefore play a minor role for the grav- 
itational potential). 

The explosion energy as a function of time is shown in 
Fig. 38. Its increase with time has several reasons: The de- 
position of energy by neutrinos is the dominant contribution 
during the early stages of the explosion after it has taken the 
LSB energy in the gain layer to a value close to zero even be- 
fore rapid expansion has begun and the explosion has set in 
(Fig. 37). When the shock reaches a radius of about 400km 
and expands farther, the post-shock matter cools sufficiently 
that neutrons and protons start recombining to ^He (or mat- 
ter accreted through the shock is not sufficiently heated to dis- 
sociate to free nucleons), which releases/saves about 7.5MeV 
per baryon and causes the energy in the gain layer to be- 
come positive. Subsequently helium recombines further to Fe- 
group elements, releasing another ~ 1 .3MeV per baryon. When 
the shock reaches about 800km (corresponding to an enclosed 
mass of 1 .45-1 .46Mq), the material accreted by the shock is no 
longer dissociated to later recombine, but instead, silicon burns 
in the shock to ^^Ni. This process, however, ends when the 
shock reaches about 1700km (1.527M0) and does not heat the 
stellar gas to Si-burning temperatures any more. At late times, 
the explosion energy is increased by the energy input from the 
baryonic wind which is driven by neutrino energy deposition 
just outside of the neutrinosphere (starting at / ^ 300ms after 
bounce). The wind material initially also fully recombines, re- 
leasing additional nuclear binding energy, but with decreasing 
densities and increasing entropy the recombination becomes 
incomplete and a large number of a-particles are left (a-rich 
freeze-out). The time evolution of the composition in the in- 
falling and exploding layers is indicated by the grey shadings 
in Fig. 36. 

On the other hand, the shock also loses energy by lifting ac- 
creted material of the outer layers in the gravitational potential 
of the star. 

The explosion energy in Fig. 38 is already positive when 
the interface between Si-layer and O-rich Si-shell passes the 
shock at 150ms. We then notice a steep rise after the onset of 
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the explosion due to neutrino heating, nucleon recombination, 
and burning. At 300ms post-bounce burning ceases and the en- 
ergy increase becomes less steep. Finally, around 400ms, the v- 
driven wind sets in, which contributes energy at a much smaller 
rate: The energy increases only very slowly. The transient little 
drop in the explosion energy between 468ms and 500ms results 
from the mapping of the model from a 2D to a ID grid. 

At the end of our simulation, the shock has reached an en- 
closed mass of \.9Mq. From Fig. 10 we can thus say that the 
shock will still lose an energy of 0.16FOE (FOE = ten to the 
fifty-one erg) for unbinding the outer layers of the star. The 
final value of the explosion energy will also be changed by fur- 
ther input from the v-driven wind, and possibly by fallback of 
material onto the PNS. Assuming that wind power and gravi- 
tational binding energy of the outer layers roughly compensate 
each other, a value of 0.5FOE seems to be a good estimate for 
the final explosion energy of our model. This is comparable 
to the energy that has been estimated for SN1999br (Hamuy 
2003) and must certainly be considered as "weak" explosion. 



The neutron star left behind by the explosion has an ini- 
tial baryonic mass of about 1.4Mo. Figure 40 shows that it is 
1.41Mo 800ms after bounce, which can decrease by some 
0.0 1M0 due to losses in the neutrino-driven wind, and later 
increase due to possible fallback. 

Mass and composition of ejecta The discussed explosion 
model reveals interesting properties also with respect to the 
neutron-to-proton ratio in the ejecta. Since we did not carry 
along marker particles in our Eulerian 2D simulation, the 
distribution of the ejecta was evaluated in a post-processing 
step: To determine the total amount of matter that is ejected 
as function of Y^, we integrate over time the flux of mass (for 
given value of Y^) through a surface of fixed radius ro, resulting 
in the expression 

M(Yi < Fe < i-z+i) = / 7^ 

COS(??min) - COS()9-max) 

n jj.cond 
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Fig. 40. Radial profiles for difi'erent quantities at post-bounce times of 200.2ms (solid), 400.6ms (dashed), 600.5ms (thick solid), and 800.9ms 
(thick dashed) for the exploding Model sl5Gio_32.a. For 200.2ms and 400.6ms laterally averaged information is shown, the two later times 
correspond to the phase of the evolution after t = 468ms which was followed in spherical symmetry. 



Here, AA = Inr^icos^j+i - cos&j) is the cell surface through 
which the matter flows and the factor in front of the sum 
renormalizes the result from the wedge used in our simulation 
to the whole star. The first sum accompUshes the integration 
over time, while the second sum integrates over those angular 
bins where the following condition is fulfilled: v,.(ro, §, t„) > 
A F; < YJro, t„) < Ym for a given bin [F,-, Y^]. 

Choosing a suitable value for ro for this evaluation is not 
a trivial task: On the one hand, it should be as large as possi- 
ble so that the matter advected through the surface at ro is no 
longer affected by neutrino processes and thus has reached 
its terminal value. On the other hand, the advection of material 
on our Eulerian grid and the associated mixing of matter tends 
to destroy small patches with extreme values of Y^. Therefore 
we would like to choose vq as smaU as possible to keep this 
numerical diffusivity as low as possible. 

The results shown in Fig. 41 were obtained with ro = 
350km. The time resolution of the processed data lies between 
one and two ms, the post-processing for Fig. 41 has ended at 



?pb = 468ms. To estimate the reliability of our analysis, we 
tested the diff'erences caused by moving ro to 300km (shift in- 
dicated by dark shading in Fig. 41) or 400km (light shading). 
In the latter case, the advection between grid zones has partly 
"washed away" small regions of, in particular, low Y^, while in 
the former case we miss some neutrino processing of matter at 
radii between 300km and 350km, which has the tendency of 
increasing Fe- The differences between the three values of ro 
can be interpreted as error bars. The integrated masses in the 
three cases witii ro = 350, 300, and 400km are Mtot = 0.0267, 
0.0273, and 0.0274Mo, respectively. 

Another problem is given by material that falls back 
through the surface at ro, and then reemerges with a difl'erent 
value of Fe; this material is counted twice in our analysis. For 
example, in Fig. 36 we can see that between 400ms and 425ms, 
some material, about IQ'^Mq, temporarily falls back below the 
radius of 350km. The amount of such double processed mate- 
rial, however, is small compared to the total integrated mass. 
Also note that the distribution of neutrino processed ejecta of 
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Fig. 41. Mass of neutiino-processed ejecta as a function of the ter- 
minal value of the electron fraction Y^. for Model sl5Gio_32.a. The 
masses are calculated with Eq. (30). The solid bars correspond to an 
evaluation of the outward directed mass flow in the 2D model through 
a sphere at a radius of tq = 350km, while the dark and light shad- 
ing indicate the variation obtained when the analysis is performed at 
ro = 300km and kq = 400km, respectively. The insert shows the mass 
distribution for values of Ye around 0.5 with increased resolution (fig- 
ure from Pruet et al. 2005). 



Fig. 41 is not the final one because it is based on an analysis 
of the explosion only until 468ms, excluding the subsequent 
neutrino-driven wind. For 7e < 0.5, however, the plot presents 
upper limits of the ejected mass at early times of the explosion, 
because the following ejecta have > 0.5 (although at very 
late times the neutrino-driven wind might become neutron-rich 
again). 

The production of proton-rich ejecta (Fe > 0.5) during the 
early phase of the supernova explosion seems to be a generic 
result of models with an elaborate spectral treatment of the 
neutrino transport either by solving the Boltzmann equation 
(Liebendorfer et al. 2003; Thielemann et al. 2003; Frohlich 
et al. 2004) or by the variable Eddington factor closure of the 
moments equations used in this work (see also Pruet et al. 
2005). According to Hoffman et al. (1996) the amount of su- 
pernova ejecta with ^ 0.47 must be less than about 10~^Mq 
per event in order to avoid overproduction of closed neutron 
shell (A^ = 50) nuclei, for example ^'^Sr, ^^Y, and ''^Zr, com- 
pared to observed galactic abundances of these elements. Our 
explosion model fulfills this constraint (Fig. 41; possible late- 
time ejection of Fg < 0.5 material in the neutrino-driven wind is 
unlikely to be sufficiently massive to be in conflict with the ob- 
servations) in contrast to previous models with partly crude or 
simple approximations of the neutrino transport (Herant et al. 
1994; Burrows et al. 1995; Janka & Mtiller 1996). 

The neutrino-processed outflow becomes proton-rich be- 
cause the neutron-to-proton mass difference favors positron 
captures, e'^-i-n — > p-i-Ve, against electron captures, e"-i-p n+ 
Ve, when the electron degeneracy is low, i.e. the electron chem- 
ical potential /j.^ ^ kT (Beloborodov 2003; Liebendorfer et al. 
2003). At large distance from the nascent neutron star, before 



weak processes freeze out, neutrinos are more energetic and 
more abundant than electrons and positrons and therefore Vg 
and Ve captures on neutrons and protons, respectively, may be 
the dominant weak reactions (in particular during the neutrino- 
driven wind phase). If the Ve and Ve spectra are not too dissimi- 
lar, the neutron-to-proton mass difference as well as weak mag- 
netism and recoil corrections of the neutrino absorption cross 
sections (Horowitz & Li 1999) favor a value of Ye ^ 0.5 (by 
reducing the Ve absorption opacity relative to the one of Ve, see 
Sect. 3.1.2, Figs. 21, 22 and Appendix A, Fig. A.l) before the 
weak processes finally freeze out when the expansion timescale 
becomes shorter than the reaction timescales. The correspond- 
ing nucleosynthesis in these neutrino-processed ejecta has been 
recently investigated in more detail for our Model sl5Gio_32.a 
by Pruet et al. (2005) and independently for ID explosion mod- 
els by Frohlich et al. (2004). 

An interesting number is the amount of ^^Ni produced by 
the explosion. Since the simulation was run with a strongly 
simplified treatment of NSE and nuclear burning, and not 
with a nuclear reaction network, we can only come up with 
a crude estimate here. Integrating the mass of all iron-group 
nuclei at the end of the simulation yields Mheavy = 0.1 2Mq. 
This value contains the ^^Ni produced by ^**Si burning in the 
shock-heated layers as well as the material ejected from the 
neutrino-processed region which recombines to iron-group nu- 
clei. Adding the 0.009Mq of ^^Ni produced by oxygen burning 
in the O-rich Si-layer we obtain a total estimate of ~ 0.1 3Mq 
of ^^Ni. Since we assume here that Si burns to Ni completely 
and that Ni is the favored nucleus produced by nucleon recom- 
bination in the ejecta (ignoring isotopes and other elements), 
our number is certainly an upper limit for the true ^^Ni yield. 
Some of this matter is also likely to fall back to the PNS instead 
of being ejected in our fairly weak explosion, reducing the final 
Ni yield even more below our estimate. 

Neutrino emission The evolution of the neutrino luminosities 

and mean energies is plotted in Fig. 42. We can distinguish 
four phases; The collapse phase, the prompt burst of electron 
neutrinos, the accretion phase and, after the onset of the explo- 
sion, the beginning of the Kelvin-Helmholtz neutrino cooling 
phase of the nascent neutron star. The latter is marked by the 
sudden and pronounced drop of the luminosities which occurs 
when the accretion of matter onto the PNS is terminated by the 
start of the explosion. During this early phase of PNS cooling 
the neutron star deleptonizes, contracts, and heats up. For this 
reason the mean energies of the emitted neutrinos are still ris- 
ing and the neutrino luminosities level off to a nearly constant 
value with no sign of long-time decay. We, however, point out 
here that our results during this phase depend on the employed 
approximation for relativistic gravity. The effective potential 
introduced by Rampp & Janka (2002) turned out to be a bit 
too strong and to overestimate the PNS contraction and heating 
and thus the energies and luminosities of the emitted neutrinos 
(see Liebendorfer et al. 2005). It also can not ensure the correct 
Newtonian limit at large radii and thus leads to an overestima- 
tion of the gravitational force and fluid acceleration in the in- 
fall region ahead of the shock. This may be the reason why the 
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Fig. 42.Evolution of the luminosities for 
Ve (solid), Ve (dastied), and v,^ (dotted) 
in Model sl5Gio_32.a around core bounce 
(left) and during the later post-bounce 
phases (right), evaluated at a radius of 
400km for an observer at rest. For the post- 
bounce phase, we also plot the average en- 
ergies of the emitted neutrinos. 
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Fig. 43.Luminosity spectra for Ve, Vc, and Vr 
for the 2D Models sl5Gio_32.a (right) and 
sl5Gio_32.b (left) for an observer at rest 
at 400km. The spectra (given as averages 
over all directions) are shown for three dif- 
ferent post-bounce times: Shortly after the 
peak of the Ve shock breakout burst, when 
the shock reaches its maximum radius in 
Model sl5Gio.32.b, and shorfly after Model 
sl5Gio_32.a has exploded. 



mass accretion rate at late post-bounce times (in nonexploding 
models) reveals a slight tendency of increase (Fig. 17). The im- 
provement of the treatment of gravity suggested by Marek et al. 
(2005) should perform better in this respect and might modify 
the results during the long-time PNS cooling to some degree. 

Also note the transient drop of the neutrino luminosity 
shortly before the Ve burst sets in (around the bounce time), 
see Fig. 19, upper panel, and Fig. 42, left panel. A detailed 
analysis of the numerical simulations reveals that this feature 
correlates with the moment of the most extreme compression 
of the homologous core and the onset of shock formation. The 
luminosity minimum is caused by a drop of the density and 
neutrino emission in the semitransparent regime exterior to the 
shock formation radius and thus in the non-homologously col- 



lapsing layers near the mass shell where the infall velocities 
reach the maximum value. 

The evolution of the neutrino spectra also changes due to 
the onset of the explosion. Before the explosion sets in, the 
spectra become gradually harder with time (Fig. 43). This is 
a consequence of the contraction of the PNS and the persist- 
ing inflow and compression-heating of fresh material accreted 
through the shock. These structural changes lead to an in- 
crease of the temperature at the radii at which the neutrinos of 
given energy decouple from the medium. The increase of this 
frequency-dependent neutrinospheric temperature is stronger 
for high-energy neutrinos, leading to an enhanced flux at high 
energies and thus to a harder spectrum. 

While in Model sl5Gio_32.b this tendency continues until 
late times, the spectral evolution changes in Model sl5Gio_32.a 
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when the explosion sets in: The mass accretion onto the PNS is 
abruptly stopped. Now that the inflow of compression-heated 
matter has ceased, the layer near the PNS surface where the 
high-energy neutrinos decouple quickly cools and the "neu- 
trinospheric" temperature for high-energy neutrinos does not 
increase as much as with ongoing accretion (it does not drop 
either because the energy and lepton number loss lead to con- 
traction). The neutrino spectra no longer harden. This is not in 
contradiction with the gradual rise of the average neutrino ener- 
gies which are determined by the spectra around their maxima 
(e ^ 20MeV). The details of this differential behaviour in dif- 
ferent spectral bands are sensitive to the stellar structure in the 
neutrino decoupling layer and to the interaction rates and there- 
fore depend also on the neutrino flavor. A deeper analysis goes 
beyond the scope of this work. 

The violent hot bubble convection behind the shock also 
affects the neutrino emission. When strong downflows hit the 
coohng region, transient lateral variation of the electron neu- 
trino and antineutrino flux of up to a factor two can occur 
(Fig. 44a), corresponding to a total luminosity increase in Vg 
or Ve by up to about 20%. Note that our "ray-by-ray" trans- 
port scheme tends to overestimate the anisotropy and lateral 
variability of the neutrino emission (flux and flux spectra, see 
Fig. 44) produced by the narrowly coUimated downflows. A 
fully two-dimensional treatment of the transport would cer- 
tainly reveal more isotropic neutrino emission even when the 
neutrinos are produced mostly in well localized regions ("hot 
spots") near the points where the downflows dive into the cool- 
ing layer (see e.g. Walder et al. 2005). 

Interestingly, sizeable lateral flux variations show up even 
long after the onset of the explosion, e.g. at 307ms, caused by 
downflows which still manage to penetrate inward towiirds the 
PNS against the expanding supernova ejecta. 

4. Summary and conclusions 

In this paper (Paper I of a series) we have described the em- 
ployed equations and numerical implementation of our "ray- 
by-ray plus" method for multi-group neutrino transport cou- 
pled to hydrodynamics in 2D models of stellar core col- 
lapse and supernova explosions. The existing code, named 
MuDBaTH, can be extended to a 3D version in a straightfor- 
ward way. Starting from the basic assumption that the neutrino 
phase space distribution is azimuthally symmetric around the 
radial direction, which implies that nonradial flux components 
are zero, we could directly build up the 2D spectral transport 
scheme from the ID code version Vertex of Rampp & Janka 
(2002). In course of this, the (9(v/c) moments equations for 
neutrino number, energy, and momentum were extended by the 
remaining terms containing lateral derivatives, which mediate 
the coupling between the difl'erent angular directions of the po- 
lar grid and depend on the lateral component of the fluid ve- 
locity. These terms are integrated in an operator splitting step. 
Closure of the set of moments equations is achieved by using 
variable Eddington factors, which are computed from the solu- 
tion of a model Boltzmann equation as in spherical symmetry 
(cf. Rampp & Janka 2002). The current code offers the options 
to do this separately for each angular direction of the grid or 
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Fig. 44. a Energy flux of Ve normalized to the average energy flux F 
as function of polar angle i? for different post-bounce times in Model 
sl5Gio_32.a. The results are evaluated for an observer at rest at a ra- 
dius of 400km, and the times are picked for showing maximal an- 
gular variations, b For t = 179ms after bounce the Ve flux spectra 
are shown for the lateral directions with maximal and minimal fluxes, 
again given for an observer at rest at a radius of 400km. Here, dL'^°/de 
is the spectrum corresponding to the "isotropic equivalent luminosity" 
V", calculated by assuming that the flux in one polar direction were 
representative for all other directions. 

- CPU-time saving - once on a laterally averaged stellar back- 
ground. 

In addition to the velocity-dependent terms with lateral 
derivatives in the transport part of the code, we also take into 
account the lateral components of neutrino pressure gradients 
in the hydrodynamics equations at optically thick conditions. 

We provided a compilation of tests for our numerical 
scheme and input physics, by which we attempted to assess 
potential limits of its apphcability or accuracy in the context 
of core collapse calculations. In particular the following points 
were addressed: 

(1) Our approximation of general relativity (i.e., replacing 
Newtonian gravity by an effective relativistic gravitational 
potential and including redshift and time dilation in the 
neutrino transport; cf. Rampp & Janka 2002) was tested 
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against fully relativistic calculations in spherical symmetry 
(Liebendorfer et al. 2005; Marek et al. 2005). Very good 
agreement in all relevant quantities was found for cases 
with moderately relativistic fluid velocities (up to ~ 10- 
20% of the speed of light) at least until several 100 ms after 
bounce. Moreover, by 2D simulations with different polar 
wedges and lateral boundary conditions we demonstrated 
empirically that linear momentum conservation, which can- 
not be strictly guaranteed with the use of our effective 
relativistic potential in multi-dimensional simulations (see 
Marek et al. 2005), is satisfactorily well fulfilled in practi- 
cal applications. 

(2) We explored the role of the lateral derivatives in the mo- 
ments equations and of the lateral components of neu- 
trino pressure gradients in the hydrodynamics equations. 
Both turned out to be necessary extensions when non- 
equilibirium transport is applied to the multidimensional 
case and artificial fiuid acceleration and buoyancy shall 
be avoided. The latter may occur in situations where neu- 
trinos contribute significantly to the gas pressure and are 
strongly coupled to stellar medium so that their advection 
with moving fluid elements cannot be ignored. Our treat- 
ment of neutrino-hydrodynamics in two dimensions pro- 
duces results in agreement with the expectations from a 
stability analysis of the stellar fluid, which we performed 
with a "Quasi-Ledoux" criterion that tries to take into ac- 
count (in a simple way) the lepton number exchange by 
neutrino diffusion between moving fluid elements and their 
surroundings. This at least shows that our scheme does not 
artificially instigate convective activity, although, of course, 
we cannot exclude that disregarding lateral components of 
the neutrino number and energy fluxes in our code underes- 
timates neutrino diffusion and thus hampers or suppresses 
the development of doubly diffusive instabilities below the 
neutrinosphere (in particular since their occurrence seems 
to be very sensitive to the efficiency of diffusion; cf. Bruenn 
et al. 2004). 

(3) Our transport treatment, favoring radial streaming, poten- 
tially overestimates the magnitude of asymmetry and angu- 
lar contrast of the radiation field of neutrinos that stream 
off from the neutrinosphere or cooling layer. When down- 
flows or convective eddies create hot spots with enhanced 
neutrino emission, material at the involved latitudes ex- 
terior to these emission regions receives more radiation, 
whereas the exposure of neighboring lateral directions is 
underestimated. Our hydrodynamical simulations revealed 
that the corresponding angular anisotropy is nonstationary 
and varies quickly in space and time on timescales of sev- 
eral milUseconds. Short, transient, local flux enhancements 
reach typically some 10% up to occasionally a factor of 
two compared to the average value. Performing a postpro- 
cessing test by making the extreme assumption that the 
neutrino radiation field is spherically symmetric (with the 
same value of the total luminosity) we found that the total 
heating rate in the gain layer is not significantly changed 
at any time. The average heating rate is changed only by 
a few percent and the heating rate per baryon even less. 
Distinguishing in the postprocessing evaluation between 



downflows and rising hot bubbles, we determined a de- 
crease or increase, respectively, of roughly 10% in the time- 
averaged net heating. We therefore conclude that underes- 
timating the spherical symmetry of the neutrino emission is 
not very likely to produce dynamical consequences for the 
supernova evolution in our simulations. 

(4) We performed a number of ID tests to assess the conse- 
quences of an error in the EoS of Lattimer & Swesty, which 
is usually used in our core collapse simulations at NSE con- 
ditions and at densities above 6 x 10^ gcm~^. This error 
leads to an underestimation of the mass fraction of a par- 
ticles and corresponding effects in the other thermodynam- 
ical quantities Uke pressure, entropy, and temperature for 
given density, internal energy density, and electron frac- 
tion. Replacing the Lattimer & Swesty EoS below a den- 
sity of 10^^ gcm~^ by different versions of a low-density 
EoS with 4 or 17 nuclear species (neutrons, protons, a's 
and one or 14 heavy nuclei, respectively) in NSE reveals 
the expected differences of the dependent thermodynamical 
variables behind the shock, but negligibly small effects in 
the post-bounce evolution of the models, e.g., nearly iden- 
tical shock trajectories. We determined two reasons for this 
insensitivity: Firstly, the postshock entropies in our simula- 
tions with relativistic gravity are so high that a particles 
play a relatively unimportant role and therefore the EoS 
differences are rather small (of order 10%). Secondly, the 
pressure profile in the essentially hydrostatic layer behind 
the standing accretion shock is tightly constrained and thus 
determined by the gravitational field of the nascent neutron 
star on the one hand, and by the mass accretion rate and 
jump conditions at the shock on the other hand. Differences 
of density and temperature in the gain layer, of course, af- 
fect the neutrino heating, which, however, is dynamically 
not very relevant in our models which are "quite far" from 
an explosion. 

(5) Finally, we investigated the implications of wiggles in the 
entropy profile (of about 10% ampUtude), which we dis- 
covered in the deceleration layer behind a standing accre- 
tion shock, when a physical EoS is employed in which the 
energy density includes contributions from the particle rest 
masses. This numerical phenomenon seems to be linked to 
the use of a Riemann solver in our hydrodynamics code and 
turned out to be sensitive to the ratio of pressure to energy 
density fed into the Riemann solver. The artificial entropy 
wiggles disappeared after we applied an algorithmic pro- 
cedure which allows us to extract the rest-mass contribu- 
tions from the energy density of the EoS so that the true 
internal energy is evolved in the hydrodynamics part of the 
code. Comparative calculations revealed that the ugly en- 
tropy fluctuations did not lead to any significant dynamical 
effects in ID as well as 2D supernova simulations. 

In the present paper our new neutrino-hydrodynamics code 
was applied to ID and 2D simulations of stellar core collapse 
and post-bounce evolution of a 15 Mq progenitor star. Our aim 
was to investigate the effects of different input physics. We also 
explored modifications in the neutrino transport, which change 
the neutrino-matter coupling in the heating and cooling layers 
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outside of the neutrinosphere. Our findings can be summarized 
as follows: 

(a) None of our simulations with the most complete implemen- 
tation of neutrino transport yields an explosion, neither in 
spherical symmetry nor in two dimensions. The 2D models 
were computed with a lateral wedge of about 90 degrees 
around the equatorial plane and periodic boundary condi- 
tions, a setup which is similar to the one used in the first 
generation of models of convectively supported, neutrino- 
driven supernova explosions (Herant et al. 1994; Burrows 
et al. 1995; Fryer 1999). Our models with the accurate spec- 
tral description of neutrino transport and neutrino-matter 
interactions are therefore unable to reproduce the success 
of these previous simulations, in which a simphfied treat- 
ment of neutrinos by grey, flux-Umited diffusion was ap- 
plied. Since the transport is the major difference between 
the older models and our current ones (of course, there 
are other differences, too, numerical as well as in the input 
physics), we interpret our negative results as a confirmation 
of concerns raised by Mezzacappa et al. (1998b), who sus- 
pected that the transport simpUfications might have favored 
the explosions found by Herant et al. (1994) and Burrows 
etal. (1995). 

(b) Using a 90° wedge with periodic boundary conditions, 
however, imposes constraints on the size of the nonradial 
structures (i.e., the wavelengths of the modes) which are 
allowed to exist. It also prevents one from studying the 
growth of low (dipole, I - I, and quadrupole, 1 = 2) modes 
in the postshock flow as seen in neutrinoless (adiabatic) 
simulations of Blondin et al. (2003) and predicted as a con- 
sequence of the vortical-acoustic cycle (Foglizzo 2002). In 
Paper II we will demonstrate by 2D simulations with a full 
180 degree grid that a thus less constrained flow can in- 
deed develop a dominant low-mode pattern behind a highly 
deformed shock front, which performs bipolar oscillations. 
This can lead not only to quantitative but even quaUtative 
differences in the post-bounce evolution of the collapsing 
stellar core, deciding about explosion or non-explosion. 

(c) The sensitivity to the treatment of the neutrino transport be- 
came evident once more from two simulations which we 
performed without the terms in the neutrino momentum 
equation, which depend on the radial velocity (the terms re- 
sponsible for neutrino advection were, however, included). 
In contrast to the corresponding simulations with the most 
complete transport implementation these two models pro- 
duced explosions, caused by roughly a factor of two less en- 
ergy loss and more neutrino heating due to 20-30% higher 
neutrino energy densities in the cooUng and heating lay- 
ers, respectively. These layers between neutrinosphere and 
stalled shock in the ID model underwent a number of non- 
Unear, radial oscillations with growing ampUtude and fi- 
nal runaway, similar to the k mechanism of pulsational in- 
stability of stellar atmospheres. The 2D model developed 
very strong postshock convection and started to explode 
when the interface between silicon shell and oxygen-rich 
silicon shell at 1.42 reached the shock about 150 ms 
after core bounce. The net explosion energy was rather 



low, roughly 0.5 X 10^' erg. The explosion left behind a 
nascent neutron with an initial baryonic mass of about 
1.41 Mq and produced less than ~ 0.13 M© of iron-group 
elements. The latter two numbers may change by the ef- 
fects of the subsequent neutrino-driven baryonic wind and 
later fallback. Most interesting of this "artificial" explo- 
sion, however, is the electron fraction in the neutrino- 
heated ejecta. Most of this material has values above 0.5 
and less than about 10""* Mq have ^ 0.47, in agreement 
with constraints derived from observed galactic chemical 
abundances (Hoffman et al. 1996). The model therefore 
does not show the enormous overproduction of A' = 50 
closed neutron shell nuclei found in previous simulations 
(see also Pruet et al. 2005). Values of the electron frac- 
tion above 0.5 are a consequence of our accurate, spectral 
description of electron neutrino and antineutrino transport, 
enhanced by the inclusion of weak magnetism corrections 
in the charged-current reactions with nucleons. These cor- 
rections reduce the absorption of Ve relative to that of Ve 
in the neutrino-heated ejecta (cf. Horowitz & Li 1999). 
The proton richness should not depend on the manipu- 
lation of the neutrino transport which caused the explo- 
sion, because electron neutrinos and antineutrinos were af- 
fected in the same way. Moreover, values of Fe above 0.5 
were also obtained in the early phase of the neutrino-driven 
wind, which set in (at ~ 300 ms p.b.) after the supernova 
explosion had been lauched and which we followed from 
~ 470 ms to more than a second after bounce in ID with 
the velocity-dependent terms in the neutrino momentum 
equation turned on again. Finally, our result is supported by 
FrohUch et al. (2004), who used spherically symmetric gen- 
eral relativistic simulations in which the numerical treat- 
ment of the transport differed from ours and in which explo- 
sions were achieved by other modifications of the transport 
than in our work. 

(d) Our ID models reveal the existence of a (Quasi-)Ledoux- 
unstable layer below the neutrinosphere. In the correspond- 
ing 2D simulations convective activity inside the nascent 
neutron star sets in about 30 ms after bounce and persists 
until the end of the simulated evolution. Since the con- 
vective layer nearly reaches with the muon and tau neu- 
trinosphere, the luminosity of these neutrinos is enhanced. 
A detailed analysis of the consequences of this convective 
activity for neutrino emission and supernova evolution wiU 
be given in Paper II. 

(e) Improving neutrino-nucleon interactions by the effects of 
nucleon thermal motions and recoil and weak-magnetism 
corrections, and including pair annihilation and scattering 
between neutrinos of different flavors and nucleon-nucleon 
bremsstrahlung, we found important differences in many 
aspects of the transport (in spectra as well as luminosities) 
for neutrinos of all flavors. This leads to quantitative differ- 
ences of the post-bounce dynamics of the collapsing stel- 
lar core and to differences in the shock evolution, without, 
however, causing a qualitative change of the outcome of the 
simulations. 
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In Paper II we shall present simulations for different pro- 
genitors between 1 1 Mq and 25 Mq, among which we will 
also consider a case with rotation. We will explore the role of 
pre-coUapse random perturbations in the core, will study the 
effects of convection below the neutrinosphere, and will inves- 
tigate the implications of hydrodynamic instabilities - convec- 
tive and vortical-acoustic - in the layer between neutrinosphere 
and stalled supernova shock by varying the angular resolution 
and the size of the lateral wedge of the polar grid. Some results 
of these models were already published in Buras et al. (2003b) 
and Janka et al. (2004, 2005b,a). 

The spectral treatment of neutrino transport and neutrino- 
matter interactions applied in these models means a new level 
of refinement and accuracy in two-dimensional hydrodynam- 
ical supernova simulations. But our transport scheme has ad- 
vantages as well as disadvantages. On the positive side is the 
fact that it was developed on grounds of a well tested code 
for spherically symmetric problems and thus guarantees to 
produce results of well-known good quality in an important 
limit. Moreover, our neutrino-hydrodynamics code possesses 
good coarse level parallelism, which allows the efficient use 
of a large number of processors. Due to the fully impUcit 
time-differencing and the corresponding matrix inversions in 
the transport part, however, systems with a low-latency, high- 
bandwidth memory are preferable. As we argued in this paper 
and presented tests for, we believe that our transport treatment 
is a reasonably good approach for problems in which local 
macroscopic anisotropies occur but which still have an overall 
spherical geometry. The applicability of our method, however, 
becomes questionable in situations with large, global deforma- 
tion, because in this case our basic assumption of (on average) 
azimuthal symmetry of the neutrino intensity is certainly inade- 
quate. On the negative side there are also substantial CPU-time 
requirements. The disregard of neutrino flux components in the 
angular directions tends to cause an overestimation of the an- 
gular asymmetry of neutrinos streaming out from the nascent 
neutron star. This may be a handicap when accurate informa- 
tion about the corresponding directional variation is required, 
e.g., for estimating pulsar kicks or gravitational wave signals 
associated with the anisotropic neutrino emission. Future tests 
will also have to clarify the question whether the disregard of 
nonradial neutrino flux components leads to a suppression of 
doubly diffusive instabiUties inside the neutron star. A compar- 
ison with other transport methods, which might employ alter- 
native, complementary approximations to deal with the high 
dimensionality of transport in multi-dimensional problems, is 
therefore highly desirable and in fact indispensable for verifi- 
cation and validation of our approach. We are therefore looking 
forward to compare our results with those obtained by other 
groups. 
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Reaction Implementation described in References 



ve* 




ve* 


Rampp & Janka (2002) 


Mezzacappa & Bruenn (1993b); Cernohorsky (1994) 


vA 




vA 


Rampp & Janka (2002) 


Horowitz (1997); Bruenn & Mezzacappa (1997) 


vN 




vN 


Sect. A.l 


Burrows & Sawyer (1998) 


Veil 




e-p 


Sect. A.1 


Burrows & Sawyer (1999) 


VeP 




e-^n 


Sect. A.l 


Burrows & Sawyer (1999) 


VeA' 




e-A 


Rampp & Janka (2002) 


Bruenn (1985); Mezzacappa & Bruerm (1993c) 


VV 




e'e"^ 


Rampp & Janka (2002) 


Bruenn (1985); Pons et al. (1998) 


vvNN 




NN 


Rampp & Janka (2002) 


Hannestad & RafFelt (1998) 








Buras et al. (2003a) 


see Buras et al. (2003a) 


(-) (-) 
l'/<,TVe 




(-) (-) 


Buras et al. (2003a) 


see Buras et al. (2003a) 



Table A.l. Overview of all neutrino-matter interactions currently implemented in our Supernova code as our set of "improved opacities" (in 
contrast to the set of opacities listed in Rampp & Janka 2002, which we call our "standard case" containing a number of approximations 
removed in the current improved treatment). For each type of interaction we refer to the relevant section in this work or point to references 
where the fundamental aspects of the calculation of the corresponding rate are summarized and details of its numerical implementation are 
given. The third column lists references where comprehensive information can be found about the physics and the approximations employed in 
the rate calculations. In the first column the symbol v represents any of the neutrinos Ve, Ve, V/j, V/j, Vj, Vj, the symbols e", e"^, n, p, and A denote 
electrons, positrons, free neutrons and protons, and heavy nuclei, respectively. The symbol N means n or p. 



Appendix A: Neutrino opacities 

The description of the neutrino opacities which we call our "standard" set closely follows Bruenn (1985) and Mezzacappa & 
Bruenn (1993a,b) with the only exception that we, in addition, take into account neutrino pair processes due to nucleon-nucleon 
bremsstrahlung. A complete list of the considered reactions, corresponding pointers to the literature and details of the numerical 
implementation into our transport code can be found in the appendix of Rampp & Janka (2002). 

The "improved description" of neutrino opacities which is used in our latest simulations comprises the reactions summarized 
in Table A.l. Specifically, the so-called iso-energetic or elastic approximation (cf. Bruenn 1985; Reddy et al. 1998), which is 
the standard simplification for calculating rates of neutral-current neutrino scatterings off free nucleons and charged-current 
absorption reactions on nucleons, has been abandoned in order to take into account energy exchange due to nucleon recoil and 
thermal motions (e.g. Schinder 1990) as well as nucleon-nucleon correlations in the dense medium (Burrows & Sawyer 1998, 
1999; Reddy et al. 1998, 1999). Modifications of the neutrino opacities due to the weak magnetism corrections associated with 
the large anomalous magnetic moments of the proton and the neutron are also accounted for (cf. Horowitz 2002). Moreover, we 
employ in the reaction cross sections a density dependent effective mass of the nucleon, which at nuclear densities is different 
from its vacuum value, and also take into account the possible quenching of the axial- vector coupling in nuclear matter (Carter & 
Prakash 2002). In the following we shall describe the numerical handling of these neutrino-nucleon interactions in some detail. 

Note that the complete list of considered interactions also includes neutrino pair production by nucleon-nucleon 
bremsstrahlung and in particular also the flavor-coupling neutrino interactions (last two lines of Table A.l), which until re- 
cently have received only little attention from core-collapse supernova modelers. The significance of these reactions as well as 
implementation details were discussed elsewhere (Buras et al. 2003a). These processes are therefore not included in the following 
discussion. 

A.1. Neutrino-nucleon interactions 

A.1 .1 . Neutrino-nucleon scattering (vN ^ vN) 

Differential rates for inelastic scattering of neutrinos off free nucleons are calculated according to Eq. (38) of Burrows & Sawyer 
(1998). For incorporating this type of neutrino-matter interactions into our transport code the formalism developed for inelastic 
scattering of neutrinos off electrons (NES) can be exploited (see Rampp & Janka 2002, Appendix A). Different from NES, 
however, the angular integration (let (o denote the cosine of the scattering angle) of the scattering kernels /?(e, e', (o) which 
yields the coefficients 0/(e, e') of the corresponding Legendre expansion (cf. Bruenn 1985, Appendix C) cannot be performed 
analytically in the case of neutrino-nucleon scattering. Moreover, the characteristic width of the kernels as a function of in- and 
outgoing neutrino energies e, e' can be small (but finite) compared to the numerical resolution of the neutrino spectrum in the 
transport scheme, in which currently only of the order of 20 energy bins with a resolution of Ae/e ^ 0.3 can be afforded. Hence, in 
order to adequately sample the scattering kernels on such a coarse energy grid it is not sufficient to simply evaluate the functions 
(f>iie, e') at each combination of energies {ej+k, f/+i)- Instead, we introduce for each energy bin [ej, e^+i] a numerical sub-grid of 
N^. neutrino energies sj < e < ey+i (and equivalentfy for the final-state energies e'), and angle cosines -1 < w < 1 to compute 
Rie, e', CO) as given by Burrows & Sawyer (1998, Eq. 38) for all such combinations of e, e', and oj. For fixed values of e and e' we 
then numericaUy integrate R{e, e', cj), appropriately weighted with the Legendre polynomials Piia)) over angles to obtain 0;(e, e'). 
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Averaging 0;(e, e') over all subgrid energies e in energy bin [ej, ej+i] and summing over all energies e' in bin [ey/ , e^v+i ] the final, 
binned Legendre coefficients (piiej^i, ef+i) are computed. These are then employed in our neutrino transport scheme. In practice 
a six-point Gauss-quadrature scheme is used here for the angular and energy integrations. In order to correctly reproduce the low- 
temperature (in the SN case therefore low-density) and low-neutrino energy limit, where only tiny energy transfers e - e' between 
neutrinos and nucleons occur and the coefficients (pi become increasingly narrowly peaked functions of e - e', we calculate 0/ 
according to the (analytically tractable) iso-energetic approximation (Bruenn 1985; Reddy et al. 1998) if the density drops below 
10^ gem"'. Note that in this case the cross section is corrected a posteriori for nucleon recoil effects (see end of Sect. A. 1.4) 
but non-vanishing energy exchange between high-energy neutrinos and nucleons are ignored. At densities below 10^ gcm~', 
however, the typical neutrino energies are moderate ((ey) < 30MeV) and the scattering rate with nucleons becomes very low. 

Given 0""'(ey, ey) for ej < ef, detailed balance arguments (see Cemohorsky 1994) are exploited in order to compute 
(p°''\ej,ej') for ej > ey and the coefficients 0™(e;, e/) = 0°'"(e/, e,), corresponding to the inverse reaction (for details, see 
Rampp & Janka 2002, Appendix A). 

Once the Legendre coefficients 0o and 0i are known, the contribution of neutrino-nucleon scattering to the colUsion integral 
of the Boltzmann equation and its angular moments is calculated in exactly the same way as described for NES in Rampp & 
Janka (2002, Appendix A). 

Note that evaluating Eq. (38) of Burrows & Sawyer (1998) requires the knowledge of nucleon-nucleon interaction potentials. 
Within the framework pursued by Burrows & Sawyer (1998) the latter can be expressed in terms of the Fermi-liquid parameters, 
which, in turn, are directly related to macroscopic observables such as the incompressibility modulus and the symmetry 
energy 5v of bulk nuclear matter (see Reddy et al. 1999). Accordingly, we adopt values for the Fermi-liquid parameters which 
are consistent with the parameters chosen for the nuclear EoS we use for our simulations (Lattimer & Swesty 1991, with 

= 180 MeV, 5v = 29.3 MeV). 

A.1 .2. Absorption of electron-flavor neutrinos by free nucleons (veii ^ e~p, VeP ^ e"^n ) 

The calculation of the inverse mean free path l//l(e) for absorption of electron (anti)neutrinos by free neutrons (protons) is 
based on Eq. (2) of Burrows & Sawyer (1999). For the numerical implementation we employ the same techniques which we 
have already described for the neutral-current reactions in Sect. A. 1.1. Difi'erent from neutrino-nucleon scattering, however, the 
outgoing lepton is a charged lepton, which, in our context is assumed to be in local thermodynamic equilibrium within the stellar 
medium. Hence, the dependence of the interaction kernels on the energy e' of the outgoing lepton does not need to be retained 
and the individual sums over energy bins [ey, e/+i] are consequently replaced by the integral over the entire spectrum of energies 
e' of the outgoing charged lepton. Given l//l(e), the absorption opacity corrected for stimulated absorption K*{e) which enters 
our neutrino transport scheme can be calculated in a straightforward way (cf. Rampp & Janka 2002, Appendix A). 



A.1 .3. Effective mass of the nucleons 

Following Burrows & Sawyer (1999) we substitute an eff'ective mass m* for the nucleon mass m, wherever the latter appears 
expUcitly in the formalism. Notably, m* is also used for inverting the relation between particle density and chemical potential 
which is used to compute the chemical potentials of protons and nucleons (displaced by the nuclear interaction potential) for 
given number densities (cf. Burrows & Sawyer 1998, Eq. 1 1). 

For the efl'ective mass we adopt the parametric density dependence (Reddy et al. 1999) 

m(p) = - —. , (A. I) 

H- or • p/Pnuc 

where m is the vacuum mass of the nucleon and the constant a is defined by the value of the effective mass at nuclear saturation 
density (pnuc/w = 0.16 fm"^) being m*(pnuc) - 0.8 m. 



A.1 .4. Weak magnetism 

Correction factors for the iso-energetic cross-sections accounting both for the effects of weak magnetism and nucleon recoil 
effects are provided by Horowitz (2002). Both effects appear at the level of C?(e/m) with e being the neutrino energy and m the 
nucleon mass. For our purposes, we have to disentangle the contributions of both effects to the final rates in order to obtain a pure 
weak magnetism correction for the formalism described above, since the latter already includes recoil and also nucleon-nucleon 
correlations, but disregards weak magnetism. This is achieved by analytically averaging the differential rates which include weak 
magnetism as well as recoil and those for recoil alone over the scattering angle and using the resulting angle-independent ratio of 
both as a weak magnetism correction factor to the rates of Burrows & Sawyer (1998, 1999). For the charged-current interactions 
the angular reduction is performed by simply integrating the differential rates over all directions of motion of the charged lepton 
in the final state. In case of neutral-current scattering a (1 - w)- weighting for the differential rate is used which is motivated by 
the definition of the transport opacity (for details, see Horowitz 2002). 
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In detail, the Legendre coefficients for neutral-current scatterings as described in Sect. A. 1.1 are modified according to 
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is the ratio of the full correction factor (Horowitz 2002, Eq. 32) 
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to the correction factor for recoil alone (Horowitz 2002, Eq. 29) 
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In Eqs. (A.3, A.4) we have used the abbreviations 

e{e) = - 



and 



^(e) s 1 + 2e . 



(A.3) 



(A.4) 



(A.5) 



Numerical values of the couphng constants Cy, Ca and F2 can be found in Horowitz (2002, Table I). Since the values of Cy, 
Ca and F2 are different for the reactions vp ^ vp and vn ^ vn, the individual correction factors for protons and neutrons 
^Y"*^'" are weighted with the corresponding number densities «„ and Mp, respectively, to yield the final correction factor ^"'^ for 
neutral-current scatterings as defined in Eq. (A.2). This factor is then applied to the combined rate of neutrino scatterings off 
neutrons and protons in our code. 

Analogously, we modify the inverse mean free path for charged-current reactions (see Sect. A. 1.2) according to l//l(e) 
f'^(e)-lM(e), where 



'^WM.Rec 
'^Rec(^) 



(A.6) 



with 
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as given by Horowitz (2002, Eq. 22), and (ibid, Eq. 19) 
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(A.7) 



(A.8) 



The ±-symbol in Eqs. (A.3, A.7) translates to a positive sign for neutrinos and a negative sign for antineutrinos, respectively. 
In our transport code, we presently do not discriminate between v^^^ and the corresponding antiparticles. For these flavors we 
therefore adopt the arithmetic mean of the correction factors of Eq. (A.3) for neutrino and antineutrino scattering which in effect 
means that the terms with the +-symbol cancel for the heavy-lepton neutrinos. 

Recall that we switch from the description of Burrows & Sawyer (1998, 1999) to the iso-energetic approximation (Bruenn 
1985) if p < 10^ g cm"3 (see Sect. A.1.1). Consequently, in order to take into account nucleon recoil and weak magnetism also 
at these conditions the correction factors defined in Eqs. (A.2, A.6) are replaced by = (nn • <^i^ Rec ' "^wm Rec^/'-'^" "p^' 



and = X' 



WM,Rec 



, respectively, whenever the density falls below 10° g cm 



A.2. Quenching of the axial-vector coupling 

In the calculations of aU neutrino-matter interactions we replace the axial- vector coupling ^a = 1 -254 by an effective, "quenched" 
value g*j^, which depends on the baryon density according to 



slip) = 8a\1 



P 



(A.9) 



4.15(p„uc +p)J 

as suggested by Carter & Prakash (2002, Eq. 13). In effect, opacities are reduced roughly by a factor ig*^/gA)^ at densities 
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A.3. Discussion 
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Fig. A.l. Inverse mean free paths as functions of neutrino energy for absorption of (left column) and Ve (right column) by free nucleons. Note 
the different scales of the ordinates. Solid lines and symbols are drawn for the charged-current reaction rates computed according to Sect. A. 1 .2. 
Thermodynamic conditions are given in the top left comers of the plots. For given density effective values for the nucleon mass m* and the 
axial-vector coupling were calculated using Eq. (A.l) and Eq. (A.9), respectively. Thin dashed-dotted lines show the weak magnetism plus 
recoil correction of Horowitz (2002) applied to the iso-energetic cross section of Bruenn (1985), which itself is drawn as a dashed line. The 
dotted lines correspond to the recoil and thermal motion approximation of Schinder (1990). In all cases the neutrino phase space was assumed 
to be empty. Lines interpolate values computed on an energy grid of 100 points which are equidistantly spaced between and 380 MeV. For 
comparison, crosses show the corresponding rates on a geometrical energy grid with 17 bins and six-point Gauss-quadratures for the angular 
and spectral integrations (relevant for the description given in Sect. A. 1.2). This is the typical spectral resolution of our dynamical supernova 
simulations. Results obtained by using 30-point Gauss-quadratures instead are displayed by open diamonds. See the end of Sect. A.3 for a 
discussion. 
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Fig. A.2. Same as Figure A. 1 , but showing the inverse mean free paths for neutral-current scattering of neutrinos (left column) and antineutrinos 
(right column) as computed according to Sect. A. 1 . 1 (solid lines, symbols) in comparison with the conventional ("standard") description (dashed 
lines) and different approximations to the complete physics (dotted and dash-dotted). In order to obtain the quantity l//l(e) the differential 
scattering rate R{e, e', oj) was integrated over the spectrum of final neutrino energies e' and all values of the cosine o) of the scattering angle. 



Main effects: In Figures A.l, A. 2 opacities for neutrino-nucleon interactions computed with the procedures described above 
are compared with the iso-energetic approximation adopted by Bruenn (1985) and also with the description of Schinder (1990). 
The latter work approximately takes into account the reaction kinematics (recoil, thermal motions and final-state blocking of the 
nucleons) but disregards weak magnetism and nucleon-nucleon correlations. According to our core-collapse simulations for the 
15 Mq star the chosen combinations of values for the density p, temperature T, and electron fraction of the stellar medium 
are characteristic for the conditions in the gain layer, where neutrino heating is strongest (top row of Figs. A.l, A.2), the region 
where most of the neutrino luminosity is produced (middle row), and the interior of the forming proto neutron star (bottom row). 

As discussed in detail by Horowitz (2002) the effects of weak magnetism and nucleon recoil counterbalance each other for Ve 
while in the case of Ve both effects add up leading to an appreciable net reduction of the standard opacities (computed according 
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Fig.A.3. Stellar profile of a spherically sym- 
metric postbounce model (Model sl5Gio_ld.b 
at fpb = 0.257 s). Density p (solid line), tem- 
perature T (dash-dotted), and electron fraction 
(dashed) as functions of stellar radius r were 
used to calculate the cross sections as functions 
of p(r) shown in Fig. A.4. 



r [cm] 



to Bruenn 1985). This can be seen in Figs. A.l, A. 2 by comparing the bold, dashed lines ("elastic approximation") with the thin, 
dotted ("recoil only") and dash-dotted ("recoil plus weak magnetism") lines. While for given neutrino energy, the relative effects 
of the weak magnetism and recoil corrections are independent of the nucleon density, a sizeable additional reduction of both 
Ve and Ve opacities due to nucleon-nucleon correlations shows up at nuclear and supranuclear densities (see Burrows & Sawyer 
1998, 1999; Reddy et al. 1998, 1999; Horowitz 2002, for a thorough discussion of the underlying physics). 

This is also illustrated by Fig. A.4 which displays the ratio of the rates computed with the procedures described above to 
the iso-energetic approximations adopted from Bruenn (1985), for a typical stellar (post-bounce) profile. Using the density as a 
coordinate (and taking T(p) and Ydp) from the stellar profile shown in Fig. A.3) we plot this ratio for the charged-current (left 
column) and neutral-current interactions (right column) as function of density for three different values of the neutrino energy. 
Up to densities of lO'^ gcm"^ one observes a nearly constant ratio of the interaction rates, which is weakly dependent on the 
thermodynamic conditions of the stellar medium but increases for higher neutrino energies. Extant variations can be attributed 
to small but finite effects of the reaction kinematics (the latter are disregarded for p < 10^ g cm"^ where the ratios deviate from 
unity because of weak magnetism and recoil corrections of the interaction cross sections in the numerators; cf. Sects. A. 1.1, 
A. 1.4). While sizeable cross section reduction due to weak magnetism and recoil of the nucleons is visible for Ve (dashed lines) 
the aforementioned counterbalancing effect leads to only small deviations from the iso-energetic approximation in case of Ve 
(solid lines). For densities between lO'^ gcm"^ and lO'"* gem"-' the prominent peak in the adopted temperature profile of the 
stellar medium (see Fig. A.3) is reflected in the density variation of the interaction rates plotted in Fig. A.4 due to the temperature 
dependence of the reaction kinematics and nucleon-nucleon correlations. Interestingly, for low neutrino energies this leads to an 
increase of the cross sections compared to the iso-energetic approximation (top and middle row of Fig. A.4). In the considered 
density and temperature range, however, neutrinos are in thermodynamic equilibrium with the stellar medium and therefore 
precise values of the cross sections are not very relevant for neutrinos of low energies. In the high-energy part of the spectrum, 
however, which carries the significant part of the neutrino flux, our description of the neutrino-matter interactions leads to an 
overall reduction of the opacities for neutrinos and antineutrinos of all flavors (bottom row of Fig. A.4) and hence facilitates 
transport of neutrinos out of the proto neutron star. 

Numerical considerations: For plotting the solid lines in Figs. A.l, A. 2 the corresponding cross sections were evaluated on 
an equidistant energy grid consisting of 100 points. Such spectral resolution, however, can currently not be afforded in our 
dynamical supernova simulations. For comparison we therefore supplement the plots with values (shown by crosses) computed 
on the comparably coarse energy grid of the simulations (geometrical grid of 17 bins). For a wide range of conditions we find good 
overall agreement with the corresponding "reference" solutions. Notable deviations show up in regions with a steep variation of 
the opacity with the neutrino energy, which can be caused, e.g., by the final state blocking of the degenerate electron gas in the 
reaction Ve + n ^ e~ -i- p (see Fig. A. 5). Obviously, this phenomenon is not specific to the adopted physical description (see, 
e.g. the open triangles in Fig. A.5) but is simply due to the fact that the rates are calculated as averages over the finite width of an 
energy bin (cf. Sects. A.1.1, A.1.2). Increasing the number of quadrature weights for the corresponding numerical integrations 
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Fig. A.4. Ratios of cross sections as computed according to Sects. A.l, A.2 to the conventional ("standard") approximation for given neutrino 
energies as functions of density p (the input quantities T(p) and Y^ip) were obtained from the stellar profile shown in Fig. A.3). In the left colurrm 
we show the ratio for charged-current absorption of Ve (solid lines) and i>e (dashed lines). The right column shows the ratio for neutral-current 
scattering of Ve, Ve. and heavy-lepton neutrinos (dash-dotted lines) off free nucleons. According to our equal treatment of heavy-lepton neutrinos 
and antineutrinos the latter quantity is computed using arithmetic averages of the particle and antiparticle cross sections (see Sect. A. 1.4). The 
neutrino phase space was assumed to be empty. 
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Fig. A.5. Inverse mean free path for the absorp- 
tion of Ve by free neutrons for the electron- 
degenerate conditions given in the last row of 
Fig. A.l. The reference curve (solid line) was 
obtained with high overall spectral resolution 
whereas crosses show averages over the bins of 
a comparably coarse energy grid which is typ- 
ically used in our dynamical simulations. For 
comparison the same averaging procedure was 
also applied to the iso-energetic approximation 
of Bruenn (1985). The results are shown as a 
dashed line with open triangles marking the cor- 
responding coarse-grid average values. 



within the spectral (and angular) bins by a factor of five each, i.e. employing 30-point Gauss-quadratures (open diamonds in 
Figs. A.l, A.2), we find very good agreement with the standard six-point results (crosses). Hence, the relevant discretization 
error of the dynamical simulations is dominated by the employed overall spectral resolution (relative width of energy bins) and 
not by the specific sub-gridding used for evaluating the neutrino opacities. 



Appendix B: Moments equations in three dimensions 

To order 0(v/c) of the fluid velocities (the so-called Newtonian approximation) the full three-dimensional moments equations 
are given by (Kaneko et al. 1984, correcting a number of misprints): 

££(i),V.H.i|,^.H)-;|-fi-.V,,:2:=c<» (B.) 
c Dt\pj c at c at oe oe 

P D (tl\ ^ ^ ^ I d J dfi Idfi die?) 5(eN) 



c Dt\p I c at c at c ot oe oe 

where fi = v/c. In general spherical coordinates, the moments are defined by 

J = ^ J J '^J J ^^^^ ' ^ " 4^ / J /irndQ ,N=-^J J /irnndQ . (B.3) 

Here, n = {nr, n», n^) = (cos 0, sin cos co, sin sin (o), ju = cos 0, and and o) are the two angles defining the direction of 
propagation. The functional dependences yS = /i(t, r, §, ip), J = J{t, r, ip,e), are suppressed in the notation. The moments 
equations written in coordinate form are the following: 
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c dt^ ' c dt c dt de c dt de c dt de 

\ d(eN,,^) ldpA ^ dieNnj,) Adp^ _ ^\ + djeN,.,^) / 1 dp, _ 
\ de \ dr j de \r d& r j de \rsin& dip r 
^ djeN^r,) ldp»\ djeNM,) A^^PA^ d(eN^,^) ( 1 dpj, P, cos i? 



de \ dr j de \r d& r j de \r sin ■& dtp rsin^ 

^ djeN^r,) IdpA ^ djeN^i!,) II dpA ^ d(eN,,,) I 1 dp, ^ p, ^ p^cos& \\ _ ^^^^ 
de \ dr I de \r d§ ) de \r sind^ dtp r rsin§ )] ' 

Here Pr{t, r, tp) = Vyjc, p^(t, r, tp) = Vfi/c, and p^{t, r, tp) = Vip/c. Note that using the continuity equation yields 

cDt\p) \cdt ^ I ^ 

Id ^ d ^p^ d d\^ll d(r%) ^ 1 g(sin?»;6^) ^ 1 dpA ^^^^ 

c dt dr r d-^ r dip) \r^ dr r sin dd^ r sin dip ) 

The different moments defined in Eqs. (B.3) are linked by the following relations, which are generally vahd: 

Prr + Pihf + Ptpifi — J , 

Nrrr + A^rtW + Nrtp, — , 

Nrr» + N§§§ + N§,, = Hfi , 

Nrrq> + N^^p + = . (B.9) 
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Also note that P and N are completely symmetric tensors, i.e. the indices can be permutated arbitrarily. The number of indepen- 
dent variables therefore is #(/, //, P, N) = (1,3, 5, 7) after applying Eqs. B.9. 

In two dimensions the symmetry relation - /(w,//) holds which leads to 



(B.IO) 



and therefore #(J,H,P,N) = (1,2,3,4); note that the last equation in (B.9) becomes redundant. In one dimension one has 
Iiio,lj) - so that in addition 

H§ = Pr& = Nrr§ = Nfffffi = Nfi^ = 0. (B.ll) 

This time, the third equation in (B.9) becomes redundant. Furthermore, one obtains two new equations from comparing the 
relations in Eqs. (B.3): 

which results in #(J,H,P,N) = (1, 1, 1, 1). Now one can introduce the definitions H = Hr,K = Prr,L = Nrrr- 

In our approximation of two-dimensional transport, i.e. assuming = Eqs. (B.9-B.12) apply, but in contrast to 

the one-dimensional case not all terms with lateral gradients or components vanish. The moments equations for the remaining 
independent variables, / and H = Hr, simplify from Eqs. (B.4) and (B.5) to 



fid ^ d Pit d . , 
[cdt dr r M 



1 dir%) 1 5(sini?;8tf)\ I d{r^H) jSrdH 



dr 



-If 



c dt 



dPr Pr 1 

dr r 2r sin § 



dr ^ c dt 

1 d(m&p») 



r 2r sin # d& 



)]} 



^ \dr r 2rsini? d^ )^ \ r "^2rsini? d& j^cdt ~ 

(Id „d P&d\^^ I d(r^Pr) 1 d{sm^fia)\ dK 3K - J ^idpA PrdK 

\cdt ^ dr r d^l \r^ dr rsini? d^ I dr r \ dr ) cdt 

d ( ri 



(B.13) 



c dt 



dPr Pr 1 

dr r 2rsini? 



a(sini?j8tf)\ 



1 



^ 1 a(sinj?j8tf) \ 
r^2rsini? d& j 



1 

'c dt 



(B.14) 



With = J/e, "Tf = H/e, TC = K/e, and £, = L/e, the moments equations describing the evolution of neutrino number read 

(Id d d\^ ^( 1 ^^-2^ 



d(r%)_^ 1 



d_ 



IdPr 



c dt 



dPr Pr 



dr rsini? d§ 
1 5(sin^\ 



a(sin »p^) \ ^ j_ djr^'H) ^ pr d9< 



dt 



dr 



Id d P» d\^, 
cdt dr r M 



r Irsin^ dd 

1 dir^Pr) 



1 



d{sin&P&) \ 



r 2rsini? d^ / 



c dt 



(B.15) 



r^ dr rsint? && 



dr r 



\drj c dt 



de\ [c dt [dr r 2rsini? d& j \ 

,^ldPr_P^ 1 d(smm) \ ^Jpr^ 1 

\dr r 2r sin # M / \ 2r si 



1 di^m&Pi) 



r 2r sin & M 
disin^P^n I dp. 



sin § dd^ 



c dt 



(B.16) 



This system of moments equations (B. 1 3-B. 16) is very similar to the Newtonian, Oiyjc) moments equations in spherical symme- 
try (see Rampp & Janka 2002, Eqs. 7,8,30,31). We have set the additional terms arising from our approximative generaUzation 
to two dimensions in boldface. 



Appendix C: Specific energy witliout offset 

Here we present our new procedure how the Eulerian equations of hydrodynamics, Eqs. (1-7), together with a general EoS 
that handles NSE as well as non-NSE conditions, can be numerically treated without producing spurious local fluctuations in 
dependent thermodynamic variables like the entropy and pressure, a kind of numerical noise which was described in Sect. 2.5. 

As discussed in that section, the wiggles observed in radial profiles are linked to situations where the value of Yg - p/iep) + 1 
drops below 4/3. The Riemann solver in our PROMETHEUS code, however, expects Tg to be thermodynamically consistent with 
the adiabatic index Fad = (din p/dlnp),, in particular F^. should fulfill F^ > 4/3, which implies that p and e should go to zero at 
the same time. This, however, cannot be guaranteed with an arbitrary normalization of the energy e. On the other hand, there is 
freedom of choosing the normahzation of the specific energy in the EoS, and also the solution of the time-evolution equations of 
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a system depends only on differences of the specific energy and not on the adopted zero point of this energy. In order to employ 
an EoS which makes use of this natural freedom of the energy normalization in our hydrodynamics code, we therefore introduce 
a shift of the reference level of the energy when the latter is transferred between hydrodynamics code and EoS. 

We find that the wiggles observed in the radial profiles of dependent thermodynamic variables disappear when the total 
specific energy used in the hydrodynamics solver PROMETHEUS is redefined compared to Sect. 2.5 as 

e = eint + ekin, (C.l) 

where e-mt is the internal specific energy and eidn = \{v^ + v'^ + v^) is the kinetic energy per unit of mass. Then the EoS index 
Te = p/Ke - ekin)pl + 1 = p/[eintp] + 1 adopts values in the physically meaningful range (Fg > 4/3), which seems necessary 
for the Riemann solver used in PROMETHEUS to produce smooth hydrodynamic solutions. However, our high-density EoS, 
provided by Lattimer & Swesty (1991), employs an energy definition which includes the nucleon rest masses and the rest masses 
of the electrons, the so-called "relativistic" specific energy, plus a constant energy ofl'set: 

erel,0 = Sint + erm + eo, (C.2) 
where the total specific rest mass energy is defined as 

erm = y • — • + 7e • — • (C.3) 

mby mby 

The specific energy off'set eo is an arbitrary constant (set equal to -930.773 IMeVx 1.602 x 10"^ mby, see Lattimer & 
Swesty 1991 for details). This energy definition was used in the equation of energy solved in our hydrodynamics code hitherto, 
i.e. e = erei.o + eun. For reasons explained below, we still wish to retain the energy definition, Eq. (C.2), when evaluating the 
EoS and at the same time use the new definition Eq. (C.l) when solving the hydrodynamics. This is achieved by the following 
scheme. 

The Euler equation for the specific energy as defined in Eq. (C.l) can be derived from the Euler equation for the specific 
energy as defined in our code hitherto, Eq. (5), by subtracting the normalized sum over the Eqs. (7) and further subtracting 
Eq. (6), i.e. "(Eq. 5) - Y^ki^q. 7)kmk(p-/mby - (Eq. 6)mec'^/mby": 



,Eq. D) - Zjk^^cq. /)kmkc /OTby - (i^q- o)meC /Wby • 
:lear source term is 



Imc > 



(C.4) 



where the nuclear source term is 



j-rf "^by Mby 

The abundance of each nuclear species k evolves according to Eq. (7): 



Gnuc = - 2, —'^ Rk - —c'Qn . (C.5) 



Q 1 ^ \ d 

-TT ipYk) + ^ - ir^pYk Vr) + — (sin §pYk vi) = Ru (C.6) 

and the electron fraction according to Eq. (6): 

- (pFe) + - ^ [r^pY, Vr) + — (sin &pY, v*) = Qn • (C.7) 

at or^ ' r sin § OXT 

Starting at time step n, we calculate a hydrodynamics time step with PROMETHEUS using the initial total specific energy at 
this time, e", as defined in Eq. (C.l). Remember that PROMETHEUS only solves the LHS of the hydrodynamics equations (in 
this implementation Eqs. 1^,C.4,C.6,C.7), i.e. the source terms and gravitational eff^ects on the RHS of the equations are ignored 
for the moment. Also note that the EoS is not solved simultaneously with the hydrodynamics so that for instance composition 
changes, even in regions with NSE, originate solely from advection, Eqs. (C.6)! 

The solution of the evolution equations, after additionally applying the gravitational efi'ects and all source terms except for 
2nuc, 2n, and Rk yields s* as well as the advected nuclear abundances F* and the electron fraction Yl. The asterisk denotes that 
the nuclear reaction effects and the neutrino source term Qn of time step n + \ have not yet been taken into account. 

Next, we transform the energy to 

e:,,,=s*-e"^.Y.n-^-c'^Yl.!^.c-^e„ (C.8) 

^Tj "'by '"by 

where e'^^^ is the specific kinetic energy after the hydrodynamics step. 

This energy transformation brings us back to the energy definition originally used in Eq. (5) so that the remaining source terms 
of that equation and in the lepton number and nuclear abundance equations can now be appUed to the hydrodynamics variables. 
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Then all thermodynamic quantities and the composition can be updated by applying the EoS and nuclear burning routines exactly 
in the same way as we did in our former version of the code. This results in the new composition F^"^' and electron fraction Y^'^^ . 
We now transform the energy back to 

nr 

= etl + < - Z r ' -^-'"-^r'-^-c'-e,, (C.9) 

^Tj '"by '"by 

where we use e"^^^ = ^rei o- simple treatment, we have implicitly taken into account the nuclear source term Qnuc since 

^-1 _ . _ y (y-i _ YD .!^.c^- (7-1 - YD ■^■c-'^ -Q^^^N" , (C.IO) 

OTby mby p 

where A?" is the size of the time step n. 

In other words, because we take into account effects which transform mass into internal energy and vice versa, i.e. composition 
changes due to nuclear transmutation and transformation of electrons to massless neutrinos, at a stage where the energy is defined 
including the rest masses of the particles which change abundance, Eq. (C.2), this energy is not affected by the transmutations and 
does not change. The conversion of energy from Cint to e,m or vice versa is not visible in their sum, - ^im + erm- Although in 
NSE it is actually not necessary to advect the nuclear composition because it is unambiguously given for local thermodynamical 
conditions (p, T, Y^) by the EoS, we follow the evolution equations (Eq. C.6) for a set of discrete nuclei which reproduce the 
baryon number and charge number of the representative heavy nucleus present at NSE in the employed high-density EoS (see 
Appendix B.l in Rampp & Janka 2002, which scheme was here extended by adding one more representative nucleus), in order 
to be able to subtract the rest-mass contribution to erei,o (c.f. Eq. C.8). 

This procedure features several advantages: First, we do not need to calculate any nuclear source terms. This includes the 
NSE regime, where the nuclear composition and temperature are well defined functions of p, Y^, and erei.o so that an iteration 
between T, eiat, and the composition becomes unnecessary. Also in a regime where nuclear burning applies, the conservation 
of gjei.o ensures strict energy conservation in accordance with the philosophy of PROMETHEUS. Second, problems which may 
arise from the fact that the nucleon masses are not equal to their masses in vacuum at high densities because of nucleon-nucleon 
interactions are avoided by our implementation. For the energy transformations, Eqs. (C.8, C.9), we can safely use the vacuum rest 
masses, because we need only ensure that grei.o is correct. Its correctness is necessary to derive correct temperatures, composition, 
and the relativistic gravitational potential (via the gravitational mass, which includes terms depending on Crei = e-mt+Crm in general 
relativity). At the high densities in the PNS, however, the use of e (Eq. C.l) without detailed knowledge of the efi'ective nucleon 
masses once more leads to the problem that Ye - + 1 may adopt values which are not perfectly consistent with the 

EoS index which relates gas pressure and internal energy density. However, this discrepancy now occurs at conditions where 
r^. » 4/3 and the relative error is much smaller than at low densities. Our test runs with the procedure described here confirm 
that the hydrodynamics results do not suff'er from spurious oscillations anywhere in the SN core. 



